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INTRODUCTION 


The  first  year  stage  of  program  of  the  investigations  includes  the  work  with 
references,  development  of  experimental  set  up,  creation  of  necessary  devices,  and 
experiments  at  Tomsk  eiirport. 

It  is  well  known  that  data  obtained  from  measurements  of  turbulence  intensity, 
wind  velocity  profiles,  and  temperature  at  the  altitudes  of  10-15  km  signify  the  zones 
of  sharp  changes  near  the  tropopause.  In  this  context  it  is  essential  to  distinguish  the 
foot  print  of  the  aircraft  from  characteristic  of  nondisturbed  turbulence.  So  the  work 
with  references  are  still  continue  in  the  following  directions: 
atmospheric  characteristic  at  altitudes  of  10  - 15  km; 
turbulence  characteristic  at  altitudes  of  10  - 15  km; 

characteristic  of  aircraft  engines,  especially  for  B-747  and  IL-86  planes,  at  altitudes 
of  10- 15  km; 

characteristic  of  aircraft  foot  prints  under  conditions  of  real  flight. 

And  direct  measurements  of  turbulence  and  refraction  are  performing  in  a  plane  foot 
print  at  the  ground  level. 

We  are  assuming  to  construct  the  model  of  atmospheric  turbulence  on  the  high 
elevated  paths  along  the  long  distance  propagation. 

KEY  WORDS:  turbulence,  wake  vortex,  aircraft,  condensation  trail,  jet,  vortex, 
dispersion  and  diffusion  regimes,  model. 


Part  I.  CONDENSATION  TRAILS 


A  brief  review  of  the  papers  referring  to  gasdynamic,  physical-chemical  and  optical  properties  of 
condensation  trails  behind  high-altitude  aircraft  is  done. 

1.1.  Introduction 

In  Europe,  the  USA  and  the  world  over  (in  Japan  for  example)  the  study  of  aviation  emissions  and 
its  impact  on  the  Atmosphere  is  of  great  interest  during  the  last  years  [1.1 -1.9].  The  problem  of 
contrail  (condensation  trail  behind  aircraft)  formation  and  evolution  is  concerned  with  studying  of 
atmosphere  pollution,  in  particular  with  the  problem  of  impact  of  aviation  emission  on  the  ozone 
layer  [1.10]. 

ICAO  (International  Civil  Aviation  Organization)  reports  (1997)  that  in  1996  world- wide  aviation 
transported  1.38x10^  passengers  using  12980  commercial  jet  aircraft  and  3180  turboprop  aircraft 
[1.4].  World-Wide  air  traffic  increased  by  5  to  6  %  annually  on  overage  from  1970  to  1993,  and  by  7 
to  8  %  in  the  years  1994  to  1996.  Global  air  traffic  consumed  aviation  fuel  at  a  rate  of  130xl0'^  - 
180x10*^  g/year  1992  to  1995.  It  is  about  5  to  6%  of  all  petrol  products.  About  65%  of  the  fuel  is 
consumed  at  cruise  at  altitudes  between  10  and  13  km.  The  largest  fraction  is  consumed  by  wide- 

o 

body  aircraft  such  as  the  B-747  on  long  distance  flights.  Most  of  the  emissions  occur  between  30 
and  55°N  over  USA,  Europe  and  the  North  Atlantic.  Globally  the  fraction  of  fuel  burnt  above  the 
tropopause  has  been  estimated  to  be  about  34%.  Over  the  North  Atlantic  the  stratospheric  fraction 
of  fuel  consumption  is  about  50%  in  the  annual  mean.  Aviation  fuel  consumption  grew  by  about  3% 
annually  over  the  last  2  decades. 

The  cruise  altitude  of  Super  Sonic  Transport  (SST)  of  second  generation  would  be  15  -18  km  for 
Mach  number  2  aircraft  (such  as  Concorde)  or  18-20  km  for  Mach  2,4  aircraft  [1.2].  This  range  of 
altitudes  is  not  far  from:  the  peak  of  the  ozone  concentration  (about  24  km);  the  altitudes  of  the 
Polar  Stratospheric  Clouds  (PSC)  with  extend  over  the  poles  (about  20  km). 
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Our  concern  in  this  report  is  with  two-dimensional  laminar  steady  numerical  solutions  of  the 
Navier-Stokes  and  continuity  equations  in  open  cavity  geometries.  Our  ultimate  objective  is 
a  nonparallel  three-dimensional  global  linear  instability  analysis  of  these  flows.  Consequently, 
emphasis  is  presently  placed  on  the  accuracy  of  the  obtained  solutions,  since  the  sensitivity  of 
instability  results  on  the  details  of  the  basic  state  is  well  known.  The  novelty  of  the  spectral 
multidomain  algorithm  constructed  for  the  basic  flow  calculations  is  that  the  iterations  for 
satisfaction  of  solution  continuity  across  domain  interfaces  and  that  for  the  nonlinearity  of 
the  governing  system  of  equations  are  combined  in  a  single  step;  order-of- magnitude  savings  in 
computing  effort  are  thus  materialised  compared  with  the  Dirichlet-Neumann  approach  in  which 
the  aforementioned  two  steps  are  nested  within  each  other. 

Results  are  presented,  validating  the  proposed  algorithm  on  both  Poisson’s  equation  and 
well-known  solutions  of  the  incompressible  equations  of  fluid  flow  motion.  Subsequently  steady 
laminar  incompressible  flow  solutions  in  the  open  cavity  are  presented.  The  modifications  ex¬ 
perienced  by  the  steady  flow  pattern  in  the  presence  of  a  sizeable  object  inside  the  cavity  is 
documented  and  conclusions  concerning  the  global  linear  instability  analysis  of  either  prob¬ 
lem  are  drawn.  Issues  concerning  compressible  flow  simulations  for  the  problem  in  question  are 
highlighted. 

Work  which  resulted  from  questions  on  the  origin  of  numerical  residuals  in  the  course  of 
steady-state  numerical  simulations,  raised  in  the  course  of  work  performed  on  the  present  and 
Contract  No.  F61775-99  WE049,  is  submitted  as  an  independent  document  with  the  request  to 
have  its  suitability  for  a  US  Patent  examined  by  the  Air  Force. 
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V.  Theofilis 

1.  Introduction 


We  report  the  results  of  our  efforts  to  obtain  spectrally  accurate  solutions  pertaining  to  steady 
laminar  flow  in  two-dimensional  open  cavities.  This  work  is  a  pivotal  step  towards  the  tar¬ 
get  of  a  global  linear  instability  analysis  of  this  flow.  In  such  an  analysis  the  entire  flowfield 
obtained  herein  must  be  used  as  the  variable  coefficients  (so-called  basic  flow)  in  the  partial- 
differential-equation-based  complex  nonsymmetric  generalised  eigenvalue  problem  in  which  the 
Navier-Stokes,  continuity  and,  in  compressible  flow,  the  energy  equations  may  be  recast.  Conse¬ 
quently,  the  quality  of  the  expected  solutions  of  the  eigenvalue  problem  is  conditioned  by  that 
of  the  steady  laminar  basic  states  and  it  is  tempting  to  devote  sufficient  computing  power  to 
ensure  the  recovery  of  high-quality  basic  flow  results.  However,  current  hardware  technology 
limits  the  maximum  number  of  discretisation  nodes  on  which  the  respective  eigenvalue  prob¬ 
lem  may  be  solved.  This,  in  turn,  results  in  the  necessity  for  basic  flow  results  to  be  obtained 
that  are  as  accurate  as  possible  on  the  limited  number  of  discretisation  points  available.  This 
latter  requirement  points  at  spectral  methods  as  the  natural  choice  for  spatial  discretisation. 
Conversely,  relaxing  the  requirement  of  modest  resolutions,  one  may  obtain  benchmark  open 
cavity  solutions  at  reasonable  amounts  of  computing  eflFort  as  one  of  the  possible  extensions  of 
the  present  work. 

In  building  up  the  individual  elements  of  the  algorithm  used  for  the  solution  of  the  prob¬ 
lem  at  hand  we  have  followed  a  modular  approach  and  present  solutions  to  the  issues  arising 
in  both  incompressible  and  compressible  flow.  After  presentation  in  §  2  of  the  spectral  mul¬ 
tidomain  algorithm  utilised,  an  iterative  multidomain  Poisson  solver  which  forms  the  core  of 
our  incompressible  multidomain  direct  numerical  simulation  (DNS)  approach  is  validated  in  §  3 
against  two  steady  laminar  flow  problems  and  a  solution  recently  presented  in  the  Literature, 
obtained  by  an  alternative  spectral  discretisation  scheme.  Our  attention  is  subsequently  focussed 
on  the  incompressible  equations  of  motion  and  the  issues  of  domain  connectivity,  boundary  con¬ 
ditions  and  resolution  necessary  for  results  of  predefined  accuracy  to  be  obtained  are  discussed 
in  §  4  (a  -c )  by  reference  to  well-known  benchmark  solutions  of  the  equations  of  motion.  Open 
cavity  incompressible  solutions  at  a  variety  of  Reynolds  number  and  aspect-ratio  parameters  are 
presented  in  §  4  (d ).  As  a  demonstrator  of  the  potential  of  the  multidomain  algorithm  proposed 
we  recover  solutions  in  a  full-bay  model  of  the  open  cavity;  the  strong  qualitative  differences  of 
the  flow  pattern  compared  with  the  equivalent-sized  empty  open  cavity  is  thus  highlighted  and 
conclusions  regarding  a  possible  subsequent  global  linear  instability  analysis  of  either  problem 
are  drawn.  Issues  of  particular  concern  in  compressible  flow  are  discussed  in  §  5.  While  we  have 
solved  the  initial-condition  problem  for  compressible  simulations,  no  compressible  open-cavity 
solutions  have  been  obtained  so  far;  these  will  be  presented  in  due  course,  after  answers  to  the 
boundary-condition  related  issues  raised  herein  have  been  addressed  in  a  satisfactory  manner. 
Concluding  remarks  and  an  outline  of  the  several  potential  extensions  of  the  present  work  are 
presented  in  §  6. 


2.  A  spectral  multidomain  approach 

(a)  Spatial  decomposition 

The  choice  of  spectral  methods  has  been  made  on  the  grounds  of  their  superiority  over  finite- 
difference  or  finite- volume  methods  in  recovering  solutions  of  optimal  accuracy  at  modest  res¬ 
olutions.  The  nature  of  the  domain  in  which  we  are  interested  renders  single-domain  spectral 
calculations  impractical.  Within  the  framework  of  spectrally-accurate  solutions  two  options  for 
the  spatial  discretisation  of  an  open  cavity  are  a  spectral  multidomain  (Demaret  and  Deville 
1991;  Quarteroni  1991)  or  a  spectral  element  (Patera  1984)  domain  decomposition.  The  first 
approach  is  well-suited  for  domains  which  may  be  decomposed  in  regular  rectangular  subdo- 
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mains;  the  second  approach  may  either  be  formulated  in  rectangular  subdomains  or  based  on 
space  triangulation  and  is  thus  more  flexible  than  spectral  multidomain  in  the  reconstruction  of 
arbitrarily-shaped  domains,  albeit  in  the  presence  of  an  efficient  triangulation  algorithm.  Given 
the  regularity  of  the  geometry  of  the  open  cavity  we  use  here  spectral  multidomain  spatial  de¬ 
composition  based  on  rectangular  subdomains.  For  consistency  in  what  follows  we  refer  to  the 
N(orth),  E(ast),  S(outh)  and  W(est)  boundaries  and  enumerate  domains  in  a  clockwise  manner. 

Rectangular  open-cavity  type  of  two-dimensional  physical  domains  may  be  reconstructed  by 
means  of  the  2^  —  1  rectangular  subdomains  whose  sides  result  from  permutations  of  open 
and  closed  boundaries  less  that  subdomain  whose  four  sides  correspond  to  solid  walls.  These 
subdomains  are  depicted  in  Figure  1  and  may  be  classified  into  four  types,  equal  to  the  number 
of  rows  in  Figure  1,  according  to  the  number  and  location  of  the  open/closed  sides  in  each 
domain.  Using  clockwise  notation  we  start  from  the  N  boundary  in  each  domain  and  denote  a 
solid  boundary  by  T’  and  an  open  boundary  by  ’O’.  The  subdomain  whose  four  sides  are  open 
is  thus  denoted  by  0000  and  we  call  it  a  type  A  domain.  The  family  of  subdomains  having  three 
open  sides  has  four  members,  1000,  0100,  0010  and  0001  and  is  denoted  as  type  B.  The  six 
subdomains  having  two  open  and  two  closed  sides  are  1100,  0110,  0011,  1001,  0101  and  1010 
and  are  denoted  as  type  C  subdomains.  Finally,  the  family  of  subdomains  having  three  closed 
sides  and  denoted  as  type  D  has  four  members,  1101,  1110,  0111  and  1011.  Using  this  notation  a 
minimum  of  four  subdomains,  one  of  type  A,  two  of  type  B  (0010)  and  one  of  type  D  (0111)  may 
be  used  to  reconstruct  a  simple  open  cavity  geometry,  while  eight  subdomains  may  minimally 
be  used  to  decompose  space  in  the  case  of  an  open  cavity  including  a  protrusion,  one  of  type  A, 
five  of  type  B  and  two  of  type  D,  as  depicted  in  Figure  2. 

(6)  Physical  boundary  conditions  for  the  incompressible  equations  of  motion 

Before  addressing  the  issue  of  connectivity  of  the  domains  exposed  in  §  2  (a)  we  discuss  briefly 
the  system  of  equations  solved  in  incompressible  flow  in  order  for  the  issue  of  physical  boundary 
conditions  to  be  addressed  and  the  distinction  to  be  made  between  physical  boundary  conditions 
and  numerical  compatibility  conditions  necessary  to  close  the  system  of  equations  in  the  present 
multidomain  framework.  If  incompressible  solutions  are  sought,  the  system  composed  of  the 
relation  between  streamfunction  and  vorticity  and  the  steady  vorticity  transport  equation 
is  solved, 


VV  +  C  =  0,  (2.1) 

»^V2C+{V'yCx-V’xCy}=0,  (2.2) 

in  which  u  denotes  the  kinematic  viscosity  of  the  fluid.  The  choice  of  the  steady,  as  opposed 
to  the  unsteady  version  of  the  equations  of  motion,  was  made  so  that  global  two-dimensional 
amplified  linear  instabilities,  potentially  developing  upon  the  sought  laminar  basic  flow  solutions, 
are  not  allowed  to  manifest  themselves  in  the  dynamics  of  a  DNS  which  aims  at  recovery  of 
steady-state  laminar  basic  states  (Theofilis  1999).  Another  point  worthy  of  mention  regarding 
the  choice  of  the  vorticity-transport  equation  as  opposed  to  the  primitive  variable  formulation 
concerns  the  savings  expected  to  be  made  when  solving  two  as  opposed  to  three  equations  using 
potentially  long  expensive  iterative  procedures  for  the  satisfaction  of  solution  continuity  across 
interfaces.  From  an  academic  point  of  view  it  is  interesting  to  investigate  spectral  multidomain 
algorithms  in  the  context  of  the  vorticity  transport  equation,  since  multidomain  solutions  are 
typically  obtained  using  primitive  variables  (Quarteroni  1991). 

The  physical  boundary  conditions  associated  with  this  system  on  solid  walls  are  no-slip  and  no¬ 
penetration.  In  terms  of  the  velocity  components  along  the  spatial  directions  x  and  y,  respectively 
denoted  by  u  and  both  boundary  conditions  concern  the  streamfunction 
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ti  =  -0^  =  0,  and  V  =  —'ipx  =  0.  (2.3) 

Substantial  effort  has  been  put  into  translating  these  physical  boundary  conditions  on  bound¬ 
ary  conditions  for  vorticity  in  the  context  of  finite-difference  calculations  (e.g.  E  and  Liu  1996). 
In  the  present  calculations  on  solid  boundaries  we  impose  both  physical  boundary  conditions 
(2.3)  on  xp  and  derive  boundary  conditions  on  ^  using  (2.1). 

Physical  boundary  conditions  are  also  given  at  open  boundaries.  In  the  case  of  the  open 
cavity  defined  in  x  E  [xq,  oToo]  x  y  ^  [yo?  oc)  one  may  distinguish  between  three  types  of  physical 
boundary  conditions,  those  imposed  at  the  inflow  boundary  x  =  rco,  boundary  conditions  on 
the  computational  boundary  y  =  yoo  which  models  the  physical  boundary  y  ^  oo  and,  finally, 
boundary  conditions  at  the  outflow  boundary  x  =  Xqo-  The  inflow  boundary  conditions  for  the 
present  problem  are  given  in  terms  of  the  streamfunction  and  the  vorticity  of  flow  on  a  flat 
plate  at  a  given  Reynolds  number  Rcxq-  At  the  N  boundary  of  all  relevant  subdomains,  y  =  yoo 
flow  is  driven  by  u  =  1  and  v  =  0  is  imposed;  Dirichlet  and  Neumann  boundary  conditions  on 
are  used  to  impose  these  physical  boundary  conditions.  At  the  outflow  boundary  x  =  Xqo 
information  is  extrapolated  from  the  interior  of  the  computational  domain  on  the  E  boundary 
of  the  subdomain  whose  E  side  is  at  x  =  Xoo-  We  will  return  to  these  points  in  in  §  4  (a-c) 
and(d). 


(c)  Connectivity  of  the  domains  and  numerical  compatibility  conditions 

We  focus  the  discussion  on  the  open  cavity  geometry  presented  in  the  upper  part  of  Figure  2, 
in  which  Roman  numbering  of  the  subdomains  is  used.  We  indicate  interfaces  by  their  position 
with  respect  to  the  domain  followed  by  the  number  of  the  domain  in  brackets.  The  interfaces 
E(ii)/W(i),  E(i)/W(iii)  and  S(i)/N(iv)  between  the  non-overlapping  domains  (ii)/(i),  (i)/(iii) 
and  (i)/(iv),  respectively,  are  artificially  generated  internal  boundaries  resulting  from  the  spatial 
decomposition.  A  straightforward  iterative  algorithm  is  used  to  solve  the  governing  equations  in 
this  geometry  as  follows. 

The  physical  boundary  conditions  on  the  N,  S  and  W  boundaries  of  domain  (ii)  are  comple¬ 
mented  by  arbitrarily  chosen  Dirichlet  data  on  E(ii),  the  governing  equations  are  solved  in  this 
subdomain  and  Neumann  data  on  the  interface  E(ii)/W(i)  are  generated.  The  same  approach  is 
followed  on  domain  (iii)  in  which  (2. 1-2.2)  are  solved  subject  to  the  physical  boundary  conditions 
on  N,  E,  and  S  of  this  domain  and  arbitrary  Dirichlet  boundary  values  on  W(iii)  in  order  for 
Neumann  boundary  data  to  be  generated  on  the  interface  E(i)/W(iii).  Analogously,  the  (no-slip) 
physical  boundary  conditions  on  E,  S,  and  W  of  subdomain  (iv)  are  complemented  by  arbitrary 
Dirichlet  data  on  N(iv),  the  governing  equations  are  solved  and  Neumann  boundary  data  on 
the  interface  S(i)/N(iv)  are  generated.  The  system  (2.1-2. 2)  may  now  be  solved  in  subdomain 
(i)  subject  to  the  physical  boundary  condition  on  N(i)  and  the  numerical  Neumann  compatibil¬ 
ity  conditions  generated  from  the  solutions  in  the  adjacent  subdomains.  The  solution  delivers 
Dirichlet  data  on  E(i),  S(i)  and  W(i),  which  are  passed  on  to  the  adjacent  subdomains  (iii),  (iv) 
and  (ii)  using  relaxation  according  to 

9(")  =  +  (1  -  (2.4) 

where  (n  —  1)  and  (n)  denote  old  and  new  iteration  levels  and  a?  is  an  0(1)  relaxation  param¬ 
eter.  The  process  is  iterated  until  convergence;  at  convergence  continuity  is  attained  on  all 
interfaces.  The  algorithm  has  been  found  to  converge  rapidly  in  both  the  form  discussed  and  also 
when  Neumann  data  are  chosen  initially  in  subdomains  (ii),  (iii)  and  (iv)  and  Dirichlet  values 
are  used  for  the  subsequent  iteration  in  subdomain  (i). 

The  proposed  algorithm  is  the  streamfunction/vorticity-transport  analogon  of  the  so-called 
Dirichlet-Neumann  iterative  scheme  presented  in  the  context  of  primitive- variables  discretisation 
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of  the  equations  of  motion  by  Quarteroni  (1991).  In  the  case  of  two  subdomains  fii  and  O2  with 
a  common  interface  F  one  solves  the  independent  problems 


vVi"^  +  Ci"^  =  0, 

in  fii 

in  Qi 

^(n)  ^  ^(n)  ^  ^(n-1) 

on  F 

= 0, 

in  0.2 

in  ^2 

[  /dnr  =  /dnr^  /dnr  on  F 

and  updates  the  Dirichlet  data  on  the  interface  according  to 


Depending  on  whether  one  solves  the  steady  version  of  the  equations  of  motion,  as  the  case  is 
here,  or  whether  one  solves  the  time-accurate  problem  the  forcing  functions  are  defined  either 
by  (2.2), 


=  -KkS  -  }  =  0,  (2.7) 

=  (2.8) 

or  by  (2, 7-2. 8)  modified  by  the  terms  arising  from  a  semi-iraplicit  time-discretisation  (Canuto 
et  al.  1987).  In  the  latter  case  the  iterative  procedure  for  the  satisfaction  of  solution  continuity 
across  domain  interfaces  is  nested  into  that  for  the  time-stepping  algorithm.  The  key  difference 
between  the  Dirichlet-Neumann  approach  for  the  steady  equations  of  motion  and  that  for  a 
time-accurate  solution  is  that  in  the  former  case  we  do  not  distinguish  between  the  iteration 
necessary  for  the  satisfaction  of  solution  continuity  across  domain  interfaces  and  that  which 
must  be  performed  for  the  nonlinearity  of  the  vorticity-transport  equation. 

Refinement  of  spatial  discretisation  at  certain  interesting  parts  of  the  flowfield  (f.e.  solid- wall 
junctions)  is  attained  in  either  of  two  ways.  First,  retaining  the  spatial  decomposition  one  refines 
resolution  in  the  subdomains  adjacent  to  the  physical  location  of  interest.  Second,  the  subdo¬ 
mains  in  which  refinement  is  required  may  be  further  decomposed  into  smaller  subdomains, 
generating  additional  interfaces  within  the  original  subdomain  and  retaining  the  principle  of 
utilising  the  physical  boundary  conditions  in  order  to  generate  numerical  compatibility  condi¬ 
tions  for  the  newly-generated  interfaces.  For  example  subdomain  (iv)  of  Figure  2  (upper)  may 
be  decomposed  into  six  subdomains,  as  shown  in  Figure  3.  The  solution  here  starts  in  a  first 
step  of  a  ’forward’  iteration  by  assuming  Dirichlet  values  on  the  NUW  and  NUE  boundaries  of 
the  type  C  subdomains  (c)  and  (e),  respectively,  solving  the  system  of  governing  equations  in  (c) 
and  (e)  and  generating  Neumann  data  for  the  type  B  subdomains  (b),  (d)  and  (f).  In  a  second 
step,  the  physical  boundary  conditions  on  the  type  B  subdomains  are  complemented  by  the 
newly  generated  information  on  their  S,  E/W  and  S  boundaries,  respectively,  and  by  arbitrary 
Dirichlet  data  values  on  N(b),  W(b),  N(d),  N(f)  and  E(f)  and  the  governing  equations  are  solved. 
In  a  final  step,  the  newly  generated  Neumann  data  on  E(a),  S(a)  and  W(a)  are  complemented 
by  arbitrary  Dirichlet  data  on  N(a)  and  (2. 1-2.2)  are  solved  in  this  domain.  This  generates  a 
new  estimate  of  the  solution  derivative  along  the  entire  N(iv)  boundary  which  may  be  used  to 
solve  for  domain  (i)  in  the  global  solution  algorithm.  This  produces  new  function  values  along 
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N(f)UN(a)UN(b)  and  the  flow  of  information  may  be  reversed  in  a  ’backward’  local  iteration 
within  the  original  domain  (iv) . 

3,  Validation  of  the  multidomain  approach  on  the  Poisson  equation 

At  the  heart  of  the  numerical  solution  for  (2. 1-2. 2)  lies  a  spectrally-accurate  Poisson  solver.  Be¬ 
fore  entering  the  details  of  the  numerical  solution  of  the  Navier-Stokes  and  continuity  equations 
we  validate  the  procedure  of  §  2  (c )  on  the  Poisson  problem 

V^^  =  F{x,y)  (3.1) 

for  three  different  cases  of  integration  domain,  forcing  function  and  associated  boundary  con¬ 
ditions. 


(a)  Poisson’s  equation  in  a  rectangular  duct 

The  first  is  calculation  of  the  steady  laminar  flow  in  a  rectangular  duct  presented  by  Tatsumi 
and  Yoshimura  (1990)  in  which  the  forcing  is  provided  by  the  constant,  pressure-gradient  related 
term  jP(x,  y)  =  —2,  The  rectangular  integration  domain  in  this  case  is  defined  by  x  E  [—A,  A]  x 
y  E  [— l,l]i  where  A  is  the  duct  aspect  ratio.  The  boundaries  in  this  example  are  taken  to 
be  solid  walls,  on  which  the  viscous  boundary  conditions  are  applied.  Space  is  discretised  by 
the  union  of  three  subdomains,  the  1011  D-type  subdomain  |x  E  [—A^xl]  x  y  £  [— 1,1]|, 

the  1010  C-type  subdomain  ja;  E  [xl,xr]  x  y  e  [“1,1]}  and  the  1110  D-type  subdomain 

jx  E  [xr^AA]  X  y  E  [—1,1]}.  The  interfaces  are  chosen  at  xl  =  — 0.5A  and  xr  =  0.25A  and 
all  subdomains  are  discretised  using  Legendre  Gauss-Lobatto  collocation  points  to  resolve  the 
X—  and  y— spatial  directions  by  Nx  and  Ny  points,  respectively.  Results  for  the  streamfunction 
at  (a;L,0),  (0,0)  and  (x/2,0)  and  its  first  derivative  at  {xr^O)  and  {xr^O)  as  function  of 
resolution  is  shown  in  Table  1  for  two  aspect  ratios,  A  =  3  and  4.  Note  that  the  magnitude 
of  '0y(a:L,O),'0a;(0, 0),7/j^(0,0)  and  'ipy{xR,0)  is  <  10”^^.  Converged  results  which  agree  in  all 
significant  digits  with  those  presented  by  Tatsumi  Sz  Yoshimura  (1990)  may  be  obtained  using 
16  collocation  points  per  spatial  direction  in  each  subdomain  and  to  =  0.5.  Our  result  for  -0  in 
the  union  of  the  three  subdomains  is  presented  in  the  form  of  eleven  contours  from  0  to  max(V^)t 
in  Figure  4,  where  the  absence  of  any  sign  of  discontinuity  at  the  interfaces  is  visible. 

The  efficiency  of  the  present  algorithm  derives  from  two  facts;  firstly,  a  spectrally-accurate 
solution  is  obtained  in  this  problem  by  solving  linear  systems  defined  by  dense  matrices  whose 
leading-dimension  is  Nx  x  Ny  as  opposed  to  3  x  Nx  x  Ny  that  would  have  been  necessary 
for  a  single-domain  calculation  at  the  same  resolution.  Second,  the  small  number  of  iterations 
necessary  for  <  10“^"^,  also  displayed  in  Table  1,  justifies  the  use  of  an 

iterative  scheme  on  the  three  small  rather  than  a  direct  scheme  on  a  single  large  problem.  In 
the  latter  problem  one  order  of  magnitude  larger  memory  and  two  orders  of  magnitude  longer 
runtime,  compared  with  those  for  solution  in  a  single  subdomain,  would  have  been  necessary 
using  direct  inversion.  Representative  results  for  the  CPU  time  necessary  for  convergence  using 
a  single  PII  300MHz  processor  are  also  shown  in  Table  1.  Shown  are  two  numbers,  one  corre¬ 
sponding  to  the  total  cost  of  the  LU-factorisation  of  the  three  implicit  matrices  corresponding 
to  the  three  subdomains  and  the  iteration  itself  and  a  second  showing  the  latter  cost  alone.  It  is 
clear  that  most  of  the  time  is  spent  in  the  one-off  factorisation,  the  results  of  which  are  stored 
and  used  during  the  iteration.  The  iterative  approach  itself  is  parallelisable,  which  may  reduce 
the  cost  of  the  iteration  process  further. 

t  Contour  level  0  is  the  wall  itself,  while  a  single  point  in  the  centre  of  the  domain  corresponds  to  the  contour  ip  —  max(V’) 
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(b)  Poisson^ s  equation  in  a  grooved- channel  geometry 

The  second  validation  case  for  the  algorithm  of  §  2  (c)  is  the  Poisson  problem-  (3.1)  subject 
to  the  same  forcing  function  and  boundary  conditions  but  solved  in  the  domain  of  the  open 
cavity.  Aside  from  testing  the  performance  of  the  algorithm  on  the  actual  domain  in  which  the 
flow  solutions  are  to  be  obtained,  the  objective  here  is  firstly  to  test  the  symmetries  present 
in  the  solution  if  a  symmetric  integration  domain  is  chosen  and  secondly  to  assess  the  level 
of  discretisation  necessary  in  order  for  solutions  to  converge  within  predefined  tolerance  in  a 
domain  containing  geometric  singularities.  Space  is  discretised  by  the  union  of  the  three  D-type 

subdomains  jlOll  :  x  E  [-1,-0. 75]  x  y  E  [0, 1]|,  joill  :  x  E  [—0.5, 0.5]  x  y  E  [— 1,0]|, 
jlllO  :  X  E  [0.5, 1]  X  y  E  [0, 1]|,  with  the  C-type  subdomain  jlOlO  :  x  E  [—0.75,  —0.5]  x  y  E 

[0, 1]|,  and  the  B-type  subdomain  jlOOO  :  x  E  [—0.5, 0.5]  x  y  E  [0, 1]|.  Of  particular  interest  is 

the  performance  of  the  algorithm  in  the  presence  of  the  discontinuities  at  (xi^yL)  =  (—0.5,0) 
and  (xn^yti)  =  (0.5,0).  Results  for  the  streamfunction  'ip  and  its  first  derivatives  along  the 
spatial  coordinates  'ipx  and  'ipy  are  shown  in  Figures  5-7;  the  symmetry  of  the  solution  about  the 
geometrical  line  of  symmetry  x  =  0  is  evident.  Results  for  ‘ipx  and  ipy  at  six  a;— locations  along 
y  =  0.5  as  a  function  of  resolution,  which  was  kept  the  same  in  the  five  subdomains,  are  shown  in 
Table  2.  The  most  striking  observation  is  that  despite  the  presence  of  the  physical  discontinuities 
in  the  integration  domain  convergence  of  the  spectral  method  is  exponential,  showing  reduction 
of  residuals  by  an  order  of  magnitude  when  doubling  the  number  of  collocation  points  in  each 
spatial  direction.  A  resolution  of  16  collocation  points  per  spatial  direction  in  each  subdomain  is 
sufficient  to  produce  results  whose  relative  discrepancy  from  the  result  at  the  highest  resolution 
is  less  than  0.05%.  On  the  other  hand,  unlike  the  case  of  the  rectangular  duct,  the  presence 
of  corner  points  in  this  problem  results  in  at  least  64  collocation  nodes  being  necessary  for 
residuals  less  than  10“^  to  be  obtained.  The  same  message  is  conveyed  by  the  small  but  nonzero 
values  of  'ipx  oii  the  line  of  symmetry  a;  =  0  where  the  latter  quantity  assumes  values  of  0(10“^) 
at  the  highest  resolution  and  also  by  the  difference  between  function  or  derivative  values  at 
a  point  and  its  mirror  image  with  respect  to  the  geometric  line  of  symmetry  x  =  0.  Such 
discrepancies  in  a  basic  flow  are  well  within  the  bounds  of  acceptability  for  accurate  linear 
instability  analysis  results  to  be  obtained.  Our  conclusion  here  is  that  the  spectral  multidomain 
algorithm  of  §  2  (a  —  c )  is  able  to  handle  the  target  grooved  channel  geometry,  albeit  at  the 
cost  of  rather  large  number  of  discretisation  points;  we  will  return  to  this  point  in  the  context  of 
simulations  of  the  Navier-Stokes  and  continuity  equations.  Given  the  well-established  property 
of  a  spectral  method  to  deliver  results  of  comparable  accuracy  with  those  of  a  second-order 
accurate  finite-difference  or  finite- volume  method  at  approximately  an  order  of  magnitude  less 
number  of  discretisation  points  per  spatial  direction,  the  present  results  may  be  used  as  guidance 
for  the  number  of  discretisation  points  necessary  for  a  second-order  accurate  method  to  deliver 
results  in  the  problem  in  question  of  accuracy  comparable  with  that  of  the  present  spectral 
scheme. 


(c)  Poisson^s  equation  in  a  backward-facing  step  geometry 

The  final  validation  case  for  the  multidomain  algorithm  before  turning  to  the  equations  of 
motion  is  that  discussed  by  Black  (1997).  The  integration  domain  here  is  the  classic  backward¬ 
facing  step  geometry  and  the  novelty  compared  with  the  previous  two  validation  cases  is  the  open 
domain  boundary  conditions.  The  integration  domain  is  decomposed  into  two  B-type  subdomains 

joOlO  :  X  E  [-TT, 0]  X  y  E  [7r/2,7r]|  and  joOOl  :  a:  E  [0, 27r]  x  y  E  [0,7r/2]|  and  an  A-type 

subdomain  joOOO  :  x  E  [0, 27r]  x  y  E  [7r/2, 7r]|.  The  forcing  function  in  (3.1)  is  a  constant  equal 
to  unity  and  the  boundary  conditions  are 
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X 

y 

ip{x,y) 

[-TT,  27r] 

TT 

0 

(3.2) 

27r 

[0,7r] 

V'l  =  0 

(3.3) 

[0, 27r] 

0 

(7r/2)V6 

(3.4) 

0 

[0,7r/2] 

(7r/2)V6 

(3.5) 

[-7r,0] 

7r/2 

(7r/2)V6 

(3.6) 

— TT 

[7r/2,  tt] 

(7r-y)"(y/3-7r/12). 

(3.7) 

Only  graphical  results  are  presented  by  Black  (1997);  nevertheless,  to  the  degree  that  results 
may  be  compared  visually,  the  agreement  seen  between  hers  and  our  result  for  '0,  to  be  found  in 
Figure  8,  is  satisfactory.  Additionally  in  this  figure  we  present  the  first  derivatives  of  the  solution 
in  none  of  which  evidence  of  discontinuity  of  the  solution  at  interfaces  can  be  found.  The  function 
and  first  derivative  values  at  the  midpoints  of  the  interfaces,  (a?,  y)  =  (tt,  0)  and  (x,  y)  =  (0, 37r/4), 
may  be  found  in  Table  3.  Overall  the  solution  converges  faster  on  the  former,  compared  with 
the  latter  point.  Both  converge  exponentially,  although  a  number  of  collocation  points  higher 
than  64  per  spatial  direction  is  found  to  be  necessary  if  convergence  beyond  the  fourth  decimal 
place  is  desired.  With  the  spectral  multidomain  approach  validated  on  the  Poisson  problem  in 
nontrivial  geometries,  we  now  turn  our  attention  to  the  solution  of  the  system  (2. 1-2. 2). 


4,  The  incompressible  equations  of  motion  in  two  spatial  dimensions 

As  it  has  been  stressed  already,  in  view  of  the  global  linear  instability  analysis  which  will 
ultimately  follow  the  present  work,  restrictions  are  placed  upon  the  maximum  total  number 
of  discretisation  points  that  can  be  used  to  solve  for  the  steady  laminar  basic  flow.  Emphasis 
is  consequently  placed  here  on  assessing  the  minimum  number  of  collocation  points  by  which 
converged  steady  laminar  solutions  may  be  obtained  within  the  framework  of  a  spectral  mul¬ 
tidomain  discretisation.  Issues  which  have  to  be  tackled  for  the  solution  of  the  system  governing 
fluid  flow  motion,  additionally  to  those  discussed  above,  are  the  following.  First,  the  extension 
of  the  iterative  algorithm  for  the  satisfaction  of  solution  continuity  across  domain  interfaces  to 
systems  of  equations  must  be  validated.  Second,  the  points  raised  regarding  the  physical  bound¬ 
ary  conditions  as  opposed  to  the  numerical  compatibility  conditions  discussed  in  §  2(6  —  c) 
deserve  elaboration  in  the  context  of  fluid  flow  simulations.  Third,  physical  boundary  condi¬ 
tions  pertinent  to  open/solid  boundaries  must  be  discussed  in  the  present  spectral  multidomain 
framework.  Finally,  the  nonlinearity  of  the  system  of  the  governing  equations  introduces  in 
principle  an  additional  iteration  besides  that  necessary  for  the  satisfaction  of  continuity  of  the 
solution  across  subdomain  boundaries.  The  relative  performance  of  our  algorithm,  in  which  this 
additional  iteration  is  circumvented,  against  an  algorithm  which  uses  an  outer  iteration  for  the 
nonlinearity  (and  potentially  the  time-advancement)  and  a  nested  inner  iteration  based  on  the 
Dirichlet-Neumann  algorithm  for  the  satisfaction  of  solution  continuity  across  interfaces  must 
be  assessed.  These  issues  are  addressed  below,  where  we  expose  the  properties  of  the  proposed 
algorithm  by  discussion  of  its  performance  on  analytically  and  numerically  well-known  fluid 
flow  examples.  Steady  laminar  flows  in  a  plane  channel,  the  lid-driven  cavity  and  the  Blasius 
boundary  layer  are  discussed  before  solutions  to  the  open  cavity  problem  are  presented. 
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(a)  Plane  Poiseuille  flow 

The  first  building  block  of  the  algorithm  which  we  wish  to  validate  is  our  approach  to  the 
boundary  conditions  imposed  at  free-boundaries.  To  this  end  a  single-domain  calculation  with 
one  free-boundary,  that  at  outflow,  is  performed.  Plane  Poiseuille  flow  (PPF)  is  an  ideal  test 
case  for  this  problem  since  its  analytical  results  serve  both  as  inflow  boundary  conditions  and  as 
test  for  the  returned  numerical  solution  in  the  entire  domain.  In  view  of  the  C)  formulation 
used  herein  one  preprocessing  step  is  necessary,  that  of  provision  of  an  analytical  functional 
form. for  -0.  We  take  a  single  domain  in  which  x  E  [—2,  3]  denotes  the  arbitrarily  chosen  extent 
of  the  downstream  spatial  direction  and  y  E  [—1, 1]  is  taken  to  be  along  the  wall- normal  spatial 
direction.  Fortuitously  in  this  flow  the  physical  domain  of  interest  and  the  standard  Chebyshev 
Gauss-Lobatto  grids  (Canuto  et  al.  1987)  coincide  and  no  mapping  of  the  computational  coor¬ 
dinate  is  necessary  in  the  wall-normal  direction.  It  is  also  worthy  of  mention  is  that  the  spatial 
homogeneity  of  this  problem  is  not  exploited  by  use  of  a  Fourier  decomposition  of  flow  quantities 
along  the  streamwise  spatial  direction  but,  rather,  a  Jacobi  polynomial  expansion  which  permits 
recovery  of  arbitrary  solutions  is  used.  It  is  left  to  the  dynamics  of  the  flow  to  select  the  correct 
(homogeneous)  solution. 

Assuming  a  downstream  pressure  gradient  px  =  —2/ Re  analytically  results  in  the  classic 
parabolic  profile  for  u{x^y)  =  1  —  y^,  while  u  =  0.  The  definition  of  u  =  'ipy  may  be  integrated 
once  to  provide 

y)  =  ^  +  y  -  (4-i) 

where  the  choice  ip{x,y  =  —  1)  =  0  has  been  made.  The  flow  vorticity  is  also  known  analyti¬ 
cally, 


C{x,y)  =  2y.  (4.2) 

The  algorithm  solves  the  Poisson  problem  (2.1)  for  given  a  delayed  C  field  subject  to  the 
boundary  conditions  of  Table  4.  At  inflow  Dirichlet  boundary  conditions  for  -0,  provided  by  eq. 
(4.1),  complement  the  boundary  condition  v  =  —'ipx  =  0  while  at  outflow  the  prescription  of 
Briley  (1971)  is  used,  which  was  presented  as  a  boundary-layer  approximation.  It  is  interesting 
to  assess  its  performance  in  the  context  of  the  present  low- Reynolds  number  Navier- Stokes 
spectral  multidomain  simulations.  On  solid  walls  the  viscous  boundary  conditions  are  imposed 
on  -0  either  directly  for  u  ['^y  —  0)  or  indirectly  for  v  by  specifying  a  constant  0  value  obtained 
by  (4.1).  With  the  solution  of  0  known  at  the  new  iteration  level  (2.2)  is  also  written  as  a  Poisson 
problem  in  which  the  nonlinear  right-hand-side  forcing  may  be  constructed  using  current-level 
information  on  0  and  ^  data  which  are  delayed  by  one  iteration;  the  Poisson  problem  for  ^  may 
then  be  solved  subject  to  the  compatibility  conditions  for  (  also  shown  in  Table  4.  Note  that  only 
one  outflow  boundary  condition  is  required  on  this  quantity,  the  other  being  imposed  on  0  at 
this  boundary  and  all  physically  available  information  being  incorporated  in  the  streamfunction 
boundary  conditions.  A  variant  of  the  present  algorithm  would  solve  a  coupled  problem  for  0  and 
(  in  which  information  on  the  forcing  term  alone  would  be  delayed  by  one  iteration.  However, 
both  the  accuracy  and  the  convergence  rate  of  results  of  single  domain  calculations  are  such 
that  the  higher  cost  of  a  coupled  solution  is  not  justified. 

Single-domain  calculation  results  converge  quickly  to  the  parabolic  profile  for  the  stream- 
wise  velocity  component  and  the  linear  dependence  of  vorticity  on  the  wall-normal  coordinate 
throughout  the  entire  domain.  Most  significantly,  this  result  has  been  found  to  be  independent 
of  the  arbitrary  extent  of  the  integration  domain  in  the  streamwise  direction  x  and  the  resolution 
of  the  wall-normal  direction  y,  provided  that  the  latter  exceeds  8  collocation  points.  The  rate 
at  which  the  analytical  results  are  obtained  is  a  function  of  the  Reynolds  number,  with  lower 
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i?e— values  resulting  in  increasing  diagonal  dominance  of  the  matrix  discretising  the  Laplacian 
operator  in  (2.2)  and  quicker  convergence  compared  with  higher  Reynolds  number  results.  The 
most  significant  finding  of  this  test-case,  though,  concerns  the  boundary  conditions  treatment. 
Use  of  equation  (2.1)  was  found  to  deliver  machine-precision  agreement  with  the  analytically 
known  inflow  boundary  condition  for  results  obtained  on  this  problem  provide  justification 
for  use  of  Briley’s  outflow  boundary  treatment  in  all  subsequent  simulations. 

The  same  problem  is  tackled  next  using  our  multidomain  approach  in  which  space  is  discretised 
by  two  1010  C-type  subdomains 

{x  €  [-3, 0]  X  2/  e  [-1, 1]}  U  {a;  €  [0, 2]  X  2/  €  [-1, 1]} . 

The  boundary  conditions  of  Table  4  are  used  at  the  N,  S  and  W  boundaries  of  the  first  and 
the  N,  E  and  S  boundary  of  the  second  subdomain,  respectively.  At  the  interface  compatibility 
conditions  ensuring  continuity  of  'ip  and  its  first  derivative  close  the  Poisson  problems  to  be  solved 
for  the  streamfunction,  while  compatibility  conditions  based  on  (2.1)  are  used  for  (  in  both 
subdomains.  Again,  the  analytical  results  are  obtained  to  within  machine-precision.  Iteration 
history  results  for  the  streamfunction,  the  streamwise  velocity  component  and  the  (constant 
throughout  the  integration  domain)  wall-normal  derivative  of  the  flow  vorticity  against  the 
time-like  variable  r,  related  with  the  number  of  iterations  performed  N  and  the  relaxation 
parameter  co  through  r  =  cjAT,  are  presented  in  Figure  9;  u  and  are  presented  in  Figure  10.  As 
has  been  mentioned,  results  were  obtained  by  iteratively  exchanging  function  and  first-derivative 
information  between  domains  at  interfaces.  The  interesting  observation  here  is  that  the  approach 
of  absorbing  the  iteration  for  the  nonlinearity  within  that  for  satisfaction  of  solution  continuity 
at  interface  boundaries  saves  an  order  of  magnitude  of  computing  effort  in  comparison  with  an 
approach  which  uses  nested  iterations  for  solution  continuity  within  those  for  the  nonlinearity 
of  the  governing  equations. 

(6)  The  lid-driven  cavity 

The  next  test  case  is  the  classic  lid-driven  cavity  flow.  This  flow  serves  as  a  testbed  for  new 
numerical  algorithms  and  has  extensively  been  simulated  in  the  past  in  its  own  right  and  on 
account  of  the  analogies  between  the  lid-driven  and  the  problem  of  interest  here,  that  of  flow 
in  an  open  cavity,  in  which  the  shear-layer  at  the  lip  of  the  open  cavity  acts  as  the  lid  driving 
flow  inside  the  cavity.  Recently,  Theofilis  (2000)  has  demonstrated  excellent  agreement  between 
the  frequency  of  the  most  unstable  global  linear  mode  calculated  by  the  novel  global  linear 
instability  analysis  technique  and  that  measured  in  a  careful  series  of  experiments  (Aidun  et 
al  1991;  Benson  and  Aidun  1992).  One  of  the  objectives  of  the  present  study  is  to  compare 
the  steady  laminar  basic  flow  patterns  in  the  lid-driven  and  the  open  cavity  in  order  to  draw 
indirect  conclusions  upon  the  expected  global  linear  instability  behaviour  of  the  open  cavity 
problem  before  embarking  upon  the  computationally  intensive  global  linear  instability  analysis 
in  the  open  cavity. 

The  spectral  multidomain  algorithm  discussed  is  applied  to  the  lid-driven  cavity  problem, 
discretising  space  using  the  union  of  a  1011  and  a  1110  D-type  domain, 

{x  G  [0, 0.5]  X  y  G  [0, 1]}  U  {x  €  [0.5, 1]  x  y  G  [0, 1]} . 

Simulations  were  performed  at  several  Reynolds  numbers  and  correspondingly  increasing  res¬ 
olutions;  the  boundary  conditions  imposed  are  shown  in  Table  5.  Excellent  agreement  with 
the  benchmark  results  of  Schreiber  and  Keller  (1983)  has  been  obtained.  Interestingly,  these 
authors  provide  vortex  core  position  and  strength  results  as  a  converging  series  of  Richardson- 
extrapolated  data.  The  maximum  relative  discrepancy  of  our  calculations  at  Re  —  4000  from 
their  results  is  <  0.3%.  Such  an  agreement  could  only  be  obtained  by  using  in  excess  of  32  spec¬ 
tral  collocation  points  per  spatial  direction  on  account  of  two  factors.  First,  the  boundary  layers 
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forming  in  the  vicinity  of  all  four  solid  surfaces  and  in  particular  in  the  salient  corners  need  to 
be  adequately  resolved.  Second,  the  singularity  of  the  boundary  conditions  at  the  lid  endpoints 
requires  placement  of  sufficient  gridpoints  in  the  vicinity  of  (x,  y)  =  (0, 1)  and  (1, 1).  Plots  of  the 
streamfunction,  vorticity  and  velocity  components  at  Re  =  100,  obtained  using  48^  Chebyshev 
collocation  points  per  subdomain  can  be  seen  in  Figure  11.  Despite  the  satisfactory  agreement 
obtained  between  benchmark  results  and  the  present  calculations,  it  has  been  observed  in  this 
nontrivial  solution  of  the  equations  of  fluid  motion  that  a  relatively  large  number  of  collocation 
points  in  each  subdomain  is  required  in  order  for  convergence  of  the  algorithm  to  be  achieved. 

(c)  The  Blasius  boundary-layer 

The  Blasius  boundary-layer  flow  is  recovered  next  by  application  of  the  spectral  multidomain 
technique  proposed.  This  flow  is  interesting  in  a  twofold  manner  in  the  context  of  the  open- 
cavity  calculations.  First,  from  a  numerical  point  of  view  the  issue  of  boundary  conditions  in  an 
open  system  in  conjunction  with  the  proposed  iterative  spectral  multidomain  algorithm  must 
be  addressed.  Second,  the  flat-plate  boundary  layer  constitutes  the  flow  upstream  of  the  open- 
cavity  solutions  to  be  discussed  shortly  and  as  such  we  wish  to  ensure  that  it  is  well  captured  in 
its  own  right.  With  hindsight,  the  results  obtained  in  this  section  have  highlighted  the  potential 
sources  of  errors  particular  to  open  systems,  introduced  in  the  solution  by  application  of  the 
multidomain  algorithm  to  this  type  of  flows. 

We  have  solved  the  equations  of  motion  on  rectangular  domains 

|x  G  [Xin’i  ^  y  ^  [0iyoo]| 

with  Xout  chosen  to  be  located  subcritically  with  respect  to  Tollmien-Schlichting  (TS)  instability 
(Schlichting  1979).  Failure  to  satisfy  this  condition  is  expected  to  result  in  amplification  of 
two-dimensional  small-amplitude  wave-like  disturbances  which,  depending  on  the  length  of  the 
integration  domain,  will  manifest  themselves  in  the  DNS  in  the  inability  to  obtain  a  converged 
steady-state  (Theofilis  1999).  Our  intention  in  terms  of  the  open-cavity  solutions  to  follow  is 
to  ensure  that  the  only  source  of  unsteadiness  originates  from  linear  instability  of  the  cavity 
itself  and  not  from  additional  instabilities  of  different  physical  origin.  To  this  end  we  chose 
Uoc  =  1,  ^in  ^  0.5,  Xout  ^  Xin-^-lOOS*^  and  u  =  1.7x  10~^  in  SI  units.  This  results  in  (dimensional) 
displacement-thicknesses  of  the  boundary  layer  5^^  ^  0.005  and  5*^*  «  0.007  at  inflow  and 
outflow,  respectively;  the  corresponding  Reynolds  numbers  are  Rein  ^  295  and  Rcout  417 
while  the  well-known  critical  Reynolds  number  for  TS-amplification  is  ReJ  «  520  (Tollmien 
1929) .  The  wall-normal  direction  was  resolved  by  mapping  the  Chebyshev  Gauss- Lobatto  points 
?7  G  [— 1, 1]  onto  the  physical  coordinate  y  using 


y  =  l 


1  -y? 
l  +  s  +  y* 


(4.3) 


Here  I  =  ydyoo  -  yo)l{yoo  -  yo  -  2yc),  s  =  2//(yoo  -yo),  lialf  of  the  total  number  of  collocation 
points  in  the  wall- normal  direction  are  placed  between  the  wall  yo  =  0  and  yc  =  (yoo  -  yo)/<i, 
the  parameter  d  takes  values  between  3  and  5;  the  outflow  boundary  is  placed  at  j/oo  ~  0'06. 

After  obtaining  converged  solutions  of  the  governing  equations  using  single-domain  computa¬ 
tions  and  as  few  as  16  collocation  points  per  spatial  direction,  the  system  (2. 1-2.2)  was  solved 
in  two  B-type  0010  subdomains, 


ja;  e  [xin,Xc]  X  y  €  [0,yoo]}  U  ja:  G  [xc,Xout]  x  y  e  [0,yoo]) 


subject  to  the  following  boundary  conditions.  At  inflow  the  Blasius  solution  at  Xin  was  calculated 
using  a  Newton-Kantorowitz  spectral  iterative  approach  (Boyd  1989)  on  the  same  wall-normal 
grid  as  the  subsequent  Navier-Stokes  calculations  were  performed,  in  order  to  avoid  the  intro- 
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duct  ion  of  interpolation-related  error  when  transferring  the  boundary- layer  data  onto  the  DNS 
grid.  The  profiles  obtained  were  then  used  as  inflow  boundary  conditions  for  -0  and  At  the 
far-field  we  imposed  'ipy  =  I  and  C  =  0.  At  the  wall  the  usual  viscous  boundary  conditions 
'iP  z=z  'ipy  =  0  were  imposed  on  ^  and  none  on  C-  At  the  outflow  boundary  Xout  the  conditions 
=  0  and  (xx  =  O5  already  validated  on  the  plane  Poiseuille  flow,  were  imposed.  Continuity 
of  the  solution  and  its  derivatives  was  ensured  at  the  interface  boundary,  the  latter  placed  at 
Xf.  =  Xin  +  {xout  —  Xin)/2.5.  In  a  manner  analogous  with  the  calculations  in  the  lid-driven  cav¬ 
ity,  the  resolution  of  the  boundary  layer  requires  an  adequate  number  of  points  resolving  the 
boundary  layer.  On  the  other  hand,  it  is  interesting  in  the  context  of  the  present  calculations 
to  obtain  solutions  within  a  predefined  low  tolerance  on  as  small  a  total  number  of  points  as 
possible;  this  requirement  was  elemental  in  the  decision  to  use  spectral  methods  for  the  dis¬ 
cretisation  of  the  governing  equations.  It  is  well-known  from  single-domain  calculations  that  at 
high  resolutions  these  methods  suff’er  from  severe  CFL-related  time-step  limitations  when  time- 
accurate  solutions  are  obtained,  and  very  small  relaxation  parameters  in  the  context  of  iterative 
solution  approaches.  As  a  matter  of  fact,  the  largest  time-step  permitted  scales  with  the  fourth 
power  of  the  smallest  grid-spacing  (Ganuto  et  al.  1987)  while  an  analogous  situation  pertains  to 
the  relaxation  parameter.  Consequently,  very  high  resolution  runs  are  quite  expensive  to  obtain 
in  general.  The  cost  is  aggravated  in  the  context  of  the  present  multi-domain  calculations  on 
account  of  the  iteration  process  necessary  for  continuity  of  the  solution  at  the  interfaces  to  be 
ensured. 

With  these  considerations  in  mind,  solutions  on  grids  comprising  12^,  16^,  32^  and  48^  col¬ 
location  points  in  each  subdomain  were  obtained.  In  Figure  12  we  present  a  local  convergence 
history  of  the  -0  and  solutions  obtained  on  the  coarsest  grid  at  the  midpoint  of  the  first  inte¬ 
gration  subdomain  against  the  time-like  variable  r;  monotonic  convergence  of  both  quantities  is 
demonstrated.  Also  shown  is  the  collapsing  of  the  calculated  solution  onto  the  Blasius  profile. 
Three  stations  are  shown,  firstly  x  =  Xin-,  where  the  Blasius  solution  is  imposed  as  an  inflow 
boundary  condition  for  the  first  domain,  secondly  the  calculated  solution  at  the  domain  inter¬ 
face  x  =  Xc  and  thirdly  the  calculated  solution  at  x  =  Xout^  It  is  worthy  of  mention  that  the 
collapsing  of  the  solutions  obtained  is  quite  adequate  even  at  this  coarse  resolution. 

However,  an  interesting  situation  arose  on  grid  refinement  of  the  coarsest  resolution  results. 
While  the  12^  calculations  converged  monotonically,  prior  to  convergence  the  16^  results  exhib¬ 
ited  an  oscillatory  behaviour  in  the  downstream  subdomain,  the  amplitude  of  which  increased 
with  the  number  of  iterations  until  divergence  of  the  simulation;  the  same  behaviour  has  also 
been  observed  as  a  result  of  too  high  a  relaxation  parameter  cj,  in  which  case  instability  appears 
without  convergence  having  been  approached.  A  typical  result  of  the  in-phase  oscillations  of 
0  and  its  derivatives  in  the  latter  case  is  shown  in  Figure  13.  Several  potential  causes  of  this 
behaviour,  which  is  absent  in  the  context  of  single-domain  calculations  using  the  same  domain 
extent  and  the  same  resolution^  were  examined  and  extensive  experimentation  was  performed 
at  a  given  extent  of  the  integration  domain  with  different  grid-refinement,  iteration  parameters 
and  outflow  boundary  conditions  alternative  to  those  used  in  the  examples  presented  so  far.  As 
a  matter  of  fact  a  modest  increase  of  the  resolution  and  decrease  of  the  associated  co  had  no 
influence  neither  on  the  pattern  seen  in  Figure  13  nor  on  its  frequency,  leading  us  to  examine  the 
possibility  of  the  instability  being  of  physical  origin.  Although  in  the  parameter  range  examined 
Tollmien-Schlichting  instability  gives  rise  to  decaying  waves  alone,  if  TS  waves  were  present  in 
the  flow  they  would  manifest  themselves  at  conditions  of  linearity,  i.e.  near  convergence  of  the 
basic  flow  to  a  steady  state,  in  line  with  the  observations  made.  On  the  other  hand,  the  absence 
of  instability  in  single-domain  calculations  points  at  the  origin  of  the  oscillations  being  of  nu¬ 
merical  nature,  associated  with  the  existence  of  an  interface  in  the  integration  domain.  Further, 
it  may  be  conjectured  that  the  practically  absent  dissipation  of  a  spectral  scheme  results  in 
catastrophic  amplification  of  the  numerically  generated  instability.  Such  a  behaviour  might  have 
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been  absent  if  a  spatial  discretisation  scheme  with  inherent  dissipation  higher  than  that  of  the 
spectral  method  were  used. 

We  first  eliminate  the  possibility  of  the  outflow  boundary  conditions  being  responsible  for 
the  phenomenon  observed.  To  this  end,  we  artificially  stabilise  the  outflow  region  by  using  a 
buffer-domain  approach,  in  which  the  convective  terms  of  the  governing  equations  are  smoothly 
reduced  to  zero.  This  technique  has  been  used  successfully  by  a  number  of  investigators  dealing 
with  spatial  direct  numerical  simulation  of  transitional  and  turbulent  flows  (Spalart  1988,  Rist 
et  al  1996).  The  function  used  to  achieve  this  purpose  in  the  context  of  the  present  simulations 
is 


A(a;)  0.5  +  tanh  ^  (4*4) 

where  5  ==  Xc  +  Oi{xout  —  ^c)  with  a  determining  the  location  s  E  [xc,  Xout]  around  which  the 
hyperbolic  tangent  function  is  centred.  The  parameter  P  determines  the  slope  in  the  hyperbolic 
tangent  function  and  ensures  that  a  given  streamwise  grid  resolves  A(x)  well.  The  governing 
equations  are  then  modified,  replacing  (2.2)  by 

+  A(a:){^,Cx  -  i’xCy]  =  0,  (4.5) 

which  delivers  unphysical  solutions  for  (^,  ^)  in  the  part  x  E  [s,  Xout]  of  the  second  subdomain 
in  favour  of  stable  solutions  upstream  of  x  =  s.  While  one  wishes  to  have  5  Xout  and  a  steep 
slope  of  the  hyperbolic  tangent  function,  there  are  limitations  to  the  extend  that  the  buffer- 
domain  approach  may  be  used  such  that  it  fulfills  its  purpose  to  damp  outflow  oscillations  and 
at  the  same  time  is  not  too  wasteful  in  terms  of  the  collocation  points  devoted  to  resolving  A 
as  opposed  to  solving  the  governing  equations.  Several  experiments  performed  using  this  tech¬ 
nique  have  improved  the  quality  of  the  transient  solutions  compared  with  the  untreated  outflow 
boundary.  However,  at  modest  resolutions  the  bufiFer  treatment  failed  to  deliver  a  satisfactory 
solution  of  the  issue  of  numerical  instability  in  the  downstream  subdomain  near  convergence  of 
the  simulation. 

We  have  established  three  possible  remedies  of  the  numerical  instability  problem.  The  first, 
least  satisfactory,  solution  is  use  of  extremely  small  u;  values.  At  the  present  conditions  and  a 
resolution  of  16^  collocation  points  in  each  of  two  subdomains  a  value  uj  ^  10“^  appears  to 
alleviate  the  problem.  While  no  discontinuity  of  the  solution  at  the  domain  interface  may  be 
found  after  a  large  number  of  iterations,  traces  of  the  instability  in  question  may  still  be  found 
in  the  solution,  as  seen  in  Figure  14.  The  solution  shown  is  converged  to  a  relative  discrepancy 
between  successive  iterates  of  <  10"^.  However,  owing  to  the  cost  of  the  iterative  approach  with 
such  a  small  value  of  uj  further  iterations  aiming  at  reducing  this  tolerance  value  to  machine- 
accuracy  levels  were  not  performed  and  it  remains  unclear  whether  the  instability  in  question 
will  ultimately  manifest  itself  in  the  solution.  The  second  means  found  to  improve  the  behaviour 
of  the  solution  is  shifting  the  inflow  location  of  the  integration  domain  x^^  upstream  of  the 
location  Xin  =  0.5.  One  may  choose  to  retain  the  outflow  location  Xout  in  its  original  position 
or  shift  it  to  x^^^  <  XquU  effectively  solving  a  different  problem.  In  the  first  case  the  instability 
appears  and  eventually  destroys  the  simulation,  while  in  the  second  case  converged  solutions 
may  be  obtained  after  a  small  number  of  iterations  which  depends  on  the  location  of  the  new 
outflow  boundary  x^^^.  Since  TS-waves  are  more  stable  in  x  E  in  x  E 

the  latter  observation  may  suggest  a  coupling  between  physical  instabilities  and  the  iterative 
approach  for  satisfaction  of  solution  continuity  at  the  interface  as  being  responsible  for  the 
observed  numerical  instability. 

Finally,  the  most  effective  way  found  to  eliminate  the  numerical  instability  of  the  type  shown 
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in  Figure  13  is  to  increase  resolution  to  the  levels  used  for  the  singular  lid-driven  cavity  problem 
or  the  Poisson  problems  in  the  geometries  with  a  corner  singularity.  Results  have  been  obtained 
discretising  space  by  means  of  two  0010  B-type  subdomains  in  each  of  which  48^  spectral  col¬ 
location  points  have  been  used.  The  Dirichlet-Neumann  algorithm  discussed  was  used  at  the 
interfaces  and  no  use  of  the  buffer-domain  technique  has  been  made  (i.e.  A(a;)  —  1  in  (4,5)).  The 
convergence  history,  obtained  with  a  modest  value  of  a;  =  10”^  may  be  found  in  Figure  15  and 
the  spatial  distribution  of  the  solutions  is  shown  in  Figure  16.  No  traces  of  numerical  instability 
may  be  found  in  neither  the  solution  for  'ip  itself  nor  in  its  streamwise  derivative  'ipx  =  at  the 
interface,  while  continuity  of  u  naturally  follows  from  that  of  'll;  and  v;  an  analogous  situation 
pertains  for  the  flow  vorticity  The  collapsing  of  the  calculated  solutions  in  both  the  16^  and 
the  48^  calculations  onto  the  Blasius  profile  is  shown  in  Figure  17. 

In  summary,  it  appears  that  a  coupling  between  physical  instabilities  and  the  iterative  algo¬ 
rithm  for  the  satisfaction  of  continuity  at  interfaces  is  responsible  for  the  observed  numerical 
instability  at  modest  resolutions.  This  phenomenon  is  absent  in  single-domain  calculations,  where 
less  collocation  points  are  in  general  necessary  in  order  to  obtain  solutions  of  analogous  quality 
and  gives  rise  to  the  investigation  of  alternative  schemes  for  the  satisfaction  of  solution  conti¬ 
nuity  at  the  interfaces,  for  example  Neumann-Neumann  or  Robin  in  combination  with  either 
Dirichlet  or  Neumann  transmission  interface  conditions  as  an  extension  of  the  present  work.  In 
the  context  of  the  present  simulations  the  most  satisfactory  remedy  was  found  to  be  an  increase 
of  the  number  of  collocation  points  to  a  level  analogous  to  that  at  which  satisfactory  solutions 
were  obtained  in  the  lid-driven  cavity,  as  far  as  the  equations  of  fluid  motion  are  concerned,  and 
to  that  at  which  converged  solutions  of  the  Poisson  equation  in  the  grooved-channel  and  the 
backward-facing  step  geometry  were  obtained. 

(d)  Open- cavity  steady  laminar  solutions 

The  experience  gained  with  the  spectral  multidomain  scheme  in  the  problems  addressed  is 
next  utilised  to  obtain  solutions  in  the  open  cavity.  The  geometry  of  the  domain  including  some 
of  the  parameters  of  the  problem  is  shown  in  Figure  18.  The  fluid  is  taken  to  be  air,  with  a 
kinematic  viscosity  of  =  1.7  x  10”^  in  dimensional  SI  units  and  flow  is  taken  to  be  driven 
by  a  constant  dimensional  free-stream  velocity  Uqo  =  1  under  zero  streamwise  pressure- gradient 
conditions.  The  control  parameter  for  the  solutions  obtained  is  the  flow  Reynolds  number,  an 
increase  of  which  may  be  interpreted  as  downstream  shifting  of  the  open  cavity  at  the  same 
free-stream  velocity  or  increase  of  the  latter  with  the  cavity  being  kept  at  a  fixed  location. 
Other  parameters  on  which  the  solution  depends  are  the  lengths  and  xr  which  determine 
the  Reynolds  number  at  the  upstream  lip  of  the  cavity  and  at  the  outflow  of  the  integration 
domain,  as  well  as  the  depth  D  and  length  L  using  which  the  two  cavity  Reynolds  numbers 
Rcl  =  uLju  and  Reo  —  uDju  may  be  defined.  In  line  with  the  argumentation  used  for  the 
Blasius  flow  with  respect  to  sub  criticality  to  TS  instability  we  ensure  that  the  Reynolds  numbers 
of  the  flow  at  the  inflow  boundary  and  at  the  upstream  lip  of  the  open  cavity  is  kept  within 
the  known  bounds.  The  choice  of  the  parameters  L  and  JA,  on  the  other  hand,  has  been  shown 
experimentally  to  be  crucial  for  the  type  of  flow  state  to  be  expected,  laminar,  transitional  or 
turbulent  (Sinha  et  al.  1982).  Accordingly,  the  values  of  both  L  and  D  as  well  as  that  of  the 
outflow  location  Xout  =  Xin-\-XL-hL’\-XR  must  be  restricted,  such  that  steady  laminar  flow  exits 
the  integration  domain.  Without  attempting  to  match  their  (different)  conditions,  we  have  used 
the  works  of  Sinha  et  al.  (1982)  and  Gatski  and  Grosch  (1984)  for  some  guidance  with  respect 
to  the  choice  of  xl,  xr^  L  and  D. 

The  spectral  multidomain  algorithm  used  in  the  open  cavity  employs  a  domain  decomposition 
analogous  to  that  of  §  3  (6)  although  here  four  subdomains  have  been  used  to  discretise  space, 
two  0010  B-type,  one  0111  D-type  and  one  A-type  subdomain  as  shown  in  Figure  2  (upper). 
The  domain  mapping  in  the  wall- normal  coordinate  in  the  B-type  subdomains  is  (4.3)  used  for 
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the  calculation  of  the  Blasius  flow.  We  used  either  our  algorithm,  which  does  not  distinguish 
between  iteration  for  the  nonlinearity  and  that  for  solution  continuity  across  interfaces  or  the 
classic  Dirichlet /Neumann  approach  which  incorporates  the  latter  within  the  former  iteration. 
Convergence  of  the  solutions  in  the  present  context  has  been  ensured  in  a  threefold  manner. 
First,  continuity  of  the  solution  across  domain  interfaces  was  enforced  in  an  inner  iteration  of 
the  Dirichlet-Neumann  algorithm  until  the  jump  across  interfaces  was  of  the  level  of  machine 
accuracy.  Second,  results  of  relative  discrepancy  below  preset  tolerances  were  obtained  as  the 
grid  was  refined;  a  typical  value  used  was  a  relative  change  of  less  than  0.1%  between  two 
successive  grid  refinements.  Third,  the  outer  iteration  for  the  nonlinearity  was  pmsued  until  the 
relative  changes  in  the  solution  were  of  the  0(10“®). 

(i)  Square  cavity  solutions 

Steady-state  solutions  have  been  obtained  in  the  open  cavity  in  the  Reynolds  number  region 
examined.  Re  G  [50,5000].  The  spatial  convergence  of  the  solutions  has  been  assessed  by  grid 
refinement  using  8^  to  32^  collocation  points  in  each  subdomain.  The  iteration  process  converged 
quickly  at  all  resolutions;  the  iteration  history  using  20^  collocation  points  in  each  subdomain  is 
shown  in  Figure  19.  The  fact  that  all  solutions  were  obtained  using  w  =  0.1  is  a  demonstration  of 
the  robustness  of  the  multidomain  algorithm  proposed,  using  which  solutions  may  be  obtained 
in  less  than  10^  iterations  at  the  low  Reynolds  numbers  and  twice  that  number  at  the  higher 
Re— values.  Global  and  local  convergence  criteria  were  used;  in  Table  6  we  present  the  maximum 
value  of  u  at  {x  =  a;inf  +  xl  +  L/2,y  =  0)  where  it  may  be  seen  that  the  relative  discrepancy 
between  the  results  at  the  two  highest  resolutions  is  less  than  0.1%.  The  spatial  distribution  of 
the  incompressible  open-cavity  results  obtained  is  shown  in  terms  of  isolines  in  Figures  20-23. 
Worthy  of  mention  in  these  figures  are  a  number  of  facts.  First,  in  all  solutions  obtained  there 
are  no  signs  of  discontinuities  at  domain  interfaces.  Second,  the  outflow  boundary  conditions 
used  for  the  Navier-Stokes  solutions  discussed  so  far  perform  equally  well  in  the  present  problem, 
despite  its  strongly  elliptic  nature.  This  is  the  result  of  the  re-establishment  of  attached  lami¬ 
nar  boundary-layer  flow  at  the  outflow  boundary;  these  boundary  conditions  may  not  perform 
equally  well  in  case  separated  flow  reaches  the  outflow  boundary  as  a  result  of  further  increase 
of  the  Reynolds  number;  this  is  one  of  the  questions  of  academic  and  practical  interest,  which 
may  be  posed  for  further  investigation  in  a  possible  extension  of  the  present  work.  Third,  the 
boundary-layer  approximation  used  to  obtain  inflow  data  becomes  progressively  more  applicable 
as  the  Reynolds  number  is  increased,  while  at  low  Reynolds  numbers  the  DNS-obtained  result 
at  the  neighbourhood  of  the  inflow  boundary  quickly  departs  from  the  boundary-layer  solution 
imposed  at  inflow.  Fourth,  an  examination  of  the  flowfields  obtained  reveals  that  an  increasingly 
large  level  of  interaction  between  flow  outside  and  that  inside  the  cavity  takes  place  as  a  result 
of  an  increase  in  the  Reynolds  number;  at  the  same  time  the  shear-layer  emanating  from  the 
upstream  lip  of  the  cavity  is  strengthened,  as  may  be  visually  appreciated  in  the  vorticity  results. 

In  the  steady  solutions  obtained  no  indication  of  the  known  from  experiment  oscillation  of  the 
flow  between  the  upstream  and  the  downstream  vertical  walls  of  the  cavity  has  been  observed. 
However,  an  interesting  situation  related  to  this  phenomenon  arose  with  a  further  increase  of 
the  Reynolds  number  to  Re  =  6000;  the  result  may  be  found  in  Figure  24.  Using  the  same 
resolutions  as  in  the  lower  Reynolds  number  simulations  we  observed  that  upwards  of  20^  col¬ 
location  points  in  each  subdomain  were  necessary  for  converged  steady-state  solutions  to  be 
obtained.  At  lower  resolutions  the  simulation  became  unstable,  indicating  that  the  number  of 
discretisation  points  is  too  low  to  resolve  the  gradients  in  the  flow.  Interestingly,  though,  an  ex¬ 
plosive  instability  of  the  algorithm  occurred  only  at  the  lowest  resolution  used,  12^,  while  when 
using  16^  collocation  points  in  each  subdomain  periodic  non-physical  solutions  were  obtained. 
The  solution  then  appeared  to  oscillate  between  the  upstream  and  downstream  vertical  walls 
of  the  cavity  while  every  signal  built  with  either  local  or  global  quantities  appeared  to  pertain 
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to  a  limiting  cycle.  In  the  absence  of  the  higher  resolution  results  one  could  misinterpret  this 
behaviour  as  the  onset  of  physical  unsteadiness  in  the  cavity.  We  ensured  that  this  observation 
is  related  to  low  resolution  and  is  not  an  artifact  of  too-high  oo  values  by  experimenting  with  the 
latter  parameter.  The  results  of  Figure  24  were  qualitatively  repeated  when  u;  values  half  and 
an  order-of-magnitude  smaller  than  that  used  in  Figure  24  were  employed.  Worse,  by  experi¬ 
menting  with  the  value  of  a;  at  a  given  resolution,  16^,  we  were  able  to  obtain  solutions  which 
(in  a  time-accurate  framework)  may  be  interpreted  as  bearing  all  the  hallmarks  of  limiting-cycle 
behaviour.  The  iteration  history  obtained  using  co  =  0.01,  in  which  the  in-phase  oscillations 
of  the  solution  components  are  clearly  visible,  is  shown  in  Figure  25.  The  altogether  plausible 
result  for  the  corresponding  velocity  components  is  shown  in  Figure  26;  this  bears  little  resem¬ 
blance  to  the  true  solution  obtained  using  upwards  of  20^  collocation  points  per  subdomain, 
also  shown  in  this  figure.  Even  more  curiously,  at  the  lower  resolution  of  8^  collocation  points  in 
each  subdomain,  the  iteration  process  converged  quickly  to  a  steady-state  solution  which  may 
be  discarded  as  unphysical  only  on  the  grounds  of  comparison  with  the  high-resolution  results. 
Summarising  the  results  of  our  experimentations  beyond  Re  =  6000  we  note  that,  aside  from 
the  well-known  need  for  grid- convergence  studies  in  order  for  a  numerically  obtained  solution 
to  be  established  as  physically  relevant,  the  new  experience  gained  during  this  part  of  the  work 
is  that  low-resolution  simulations  may  produce  results  which  can  be  erroneously  interpreted  as 
being  related  to  the  known  from  experiment  physical  behaviour  of  the  open  cavity  solution.  In 
view  of  the  different  resolution  capacity  of  spectral  methods  on  one  hand  and  finite-difference 
or  finite-volume  schemes  on  the  other,  the  present  observations  may  be  relevant  to  simulations 
using  the  latter  numerical  methods  and  resolutions  substantially  higher  than  those  presently 
employed.  Caution  is  warranted  in  order  to  discern  between  physical  behaviour  and  numerical 
artifacts  in  open  cavity  flow  simulations. 

Increasing  the  Reynolds  number  further  to  Re  =  7000  we  were  unable  to  obtain  a  converged 
steady-state  solution  when  employing  up  to  32^  collocation  points  in  each  subdomain,  as  shown 
in  Figure  27.  The  iteration  history  signal  appears  to  be  nonlinear  throughout  the  simulation.  One 
of  the  iterates  of  the  solution  shortly  before  numerical  instability  destroys  the  simulation  is  shown 
in  Figure  28.  The  results  for  the  streamfunction  and  its  derivatives  seem  plausible,  although  one 
notices  the  large  extent  of  interaction  between  fluid  outside  and  inside  the  cavity,  as  well  as  the 
strong  gradient  forming  in  the  neighbourhood  of  the  downstream  corner  of  the  cavity;  the  latter 
is  an  indication  that  further  runs  to  must  be  performed  at  higher  resolutions.  Although  it  would 
be  tempting  to  proclaim  that  we  have  identified  the  critical  Reynolds  number  for  amplification 
of  two-dimensional  (/3  =  0)  global  linear  instabilities  (Theofilis  1999)  the  results  at  Re  =  6000 
warrant  caution  and  such  a  statement  may  only  be  made  after  simulations  at  higher  resolutions 
(and  associated  smaller  values  of  a;,  i.e.  simulations  which  are  altogether  substantially  more 
expensive  than  the  one  using  32^  collocation  points  in  each  subdomain)  have  been  performed. 
This  forms  one  of  the  possible  extensions  of  the  present  work,  namely  the  identification  through 
extensive  grid-refinement  study  of  the  critical  Reynolds  number  Recr,2D  for  amplification  of  two- 
dimensional  global  linear  instabilities.  A  further  avenue  which  may  be  pursued  in  continuation  of 
the  present  work  is  simulations  of  time-periodic  states  which  are  expected  to  set  in  past  Recr,2D 
as  a  consequence  of  neutrally-stable  (^  =  0)  global  linear  eigenmodes  which  interact  nonlinearly 
with  the  steady  laminar  basic  state.  The  objective  of  pursuing  simulations  at  higher  Reynolds 
numbers  is  to  broaden  the  scope  of  the  present  work  towards  Reynolds  numbers  relevant  to 
Air  Force  needs.  Subsequent  global  linear  instability  analyses  of  both  steady  and  time-periodic 
states  will  identify  two  different  mechanisms  both  of  which  are  of  relevance  to  high-speed  flow. 
From  a  technical  point  of  view,  the  extension  of  the  present  work  into  recovery  of  time-periodic 
basic  states  may  be  pursued  using  a  straightforward  time-accurate  extension  of  the  algorithms 
validated  herein. 
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(ii)  Rectangular  cavity  solutions 

Next,  we  have  turned  our  attention  to  open  cavities  in  which  L  ^  D.  We  have  considered  a 
constant  free-stream  velocity  u  =  l  and  kept  the  locations  of  the  inflow  domain  and  the  upstream 
lip  of  the  cavity  at  the  flxed  positions  Xin  =  0.5  and  Xin  +  rrjr,  =  0.5  +  205*,  respectively,  where 
5*  is  the  flat-plate  boundary  layer  thickness  at  inflow;  the  outflow  boundary  is  placed  at  505* 
past  the  downstream  vertical  wall  of  the  cavity.  Taking  as  reference  the  square  cavity  solution 
L  =  =  105*  wide  cavity  solutions  have  been  obtained  by  taking  D  =  105*  and  L/D  =  2(1)10 

while  deep  cavity  solutions  have  been  obtained  using  L  =  105*  and  DjL  =  2,3  and  4.  In 
all  simulations  space  was  discretised  using  the  same  conformable  subdomains  as  in  the  square 
cavity  simulations,  in  which  20  collocation  points  resolved  each  spatial  direction.  An  cj  =  0.1 
has  sufficed  for  converged  solutions  to  be  obtained  after  a  few  Ixundreds  of  iterations  at  all 
parameters. 

Figure  29  shows  the  dependence  of  the  maximum  of  the  streamwise  velocity  component  along 
the  interface  of  the  domains  II  and  IV.  This  quantity  reaches  its  minimum  in  the  square  cavity 
and  depends  in  a  nearly  linear  manner  on  the  aspect  ratio  L/D  as  the  latter  increases.  The 
physical  explanation  of  this  result  is  rather  straightforward,  pointing  to  the  fact  that  increas¬ 
ingly  more  fluid  from  the  free  stream  entrains  into  the  cavity  L/D  increases.  This  may  be 
appreciated  in  the  results  of  Figures  30-33,  where  the  solutions  for  'ip,  u,  v  and  ^  at  L/J9  =  3(2)9 
are  plotted.  Our  intention  here  is  to  highlight  the  qualitative  features  of  the  solutions  obtained 
and  hence  we  present  a  large  number  of  equally  spaced  contours  between  zero  and  the  maxima 
of  each  quantity  as  well  as  the  same  number  of  contours  between  the  minima  and  zero  of  the 
respective  quantities.  The  same  results  are  plotted  in  Figure  34  on  the  same  scale,  so  that  the 
quantitative  changes  resulting  from  a  change  of  the  aspect  ratio  may  be  appreciated.  At  about 
L/D  =  b  the  single  recirculation  region  inside  the  open  cavity  is  divided  in  two  large  areas  of 
recirculating  flow  attached  to  either  of  the  upstream  and  downstream  vertical  cavity  walls;  as 
the  aspect  ratio  of  the  cavity  increases  further  the  flow  pattern  in  the  neighbourhood  of  the 
upstream  and  downstream  vertical  walls  of  the  open  cavity  increasingly  resembles  that  of  the 
union  of  two  other  interesting  prototype  flows,  namely  the  backward-  and  forward- facing  step 
flow.  The  latter  may  be  better  visualised  in  Figure  35  in  which  isolines  of  the  steady  laminar 
solution  at  L/D  =  10  are  presented.  In  order  for  the  recirculation  regions  to  be  highlighted 
ten  isolines  of  positive  values  between  zero  and  the  respective  maxima  as  well  as  ten  isolines 
of  negative  values  between  the  respective  minima  and  zero  and  the  associated  colourbars  were 
generated  for  the  respective  flow  quantities.  At  these  parameters  the  attached  laminar  flow  at 
the  inflow  boundary  separates  past  the  upstream  lip  of  the  cavity  and  reattaches  at  the  lower 
wall  of  the  cavity  at  rr  Xin-hxi  H-  1.25i3  before  separating  again  shortly  before  the  downstream 
vertical  wall  of  the  cavity  is  reached.  Thereafter  a  laminar  attached  boundary  layer  develops  up 
to  the  outflow  boundary.  Comparisons  of  the  results  obtained  with  established  backward-  and 
forward-facing  step  benchmark  simulations  is  beyond  the  scope  of  the  present  investigation  and 
will  be  reported  separately  in  due  course  as  yet  another  possible  extension  of  the  present  work. 
Prom  a  numerical  point  of  view,  in  line  with  the  results  of  the  Poisson  equation  obtained  in 
§  3  (c ) ,  it  suffices  here  to  state  that  the  present  wide  cavity  steady-state  results  demonstrate 
that  the  spectral  multidomain  algorithm  presented  herein  is  capable  of  solving  for  flow  in  both 
the  backward  and  the  forward-facing  step  geometry  in  a  straightforward  manner. 

Finally,  we  turn  our  attention  to  deep  cavities  and  present  in  Figures  36-39  and  40,  respectively, 
the  qualitative  and  quantitative  features  of  steady  and  (  solutions  obtained  at  DjL  = 

1(1)4.  These  may  be  discussed  along  both  numerical  and  physical  lines.  From  a  numerical  point 
of  view,  as  in  all  open  cavity  solutions  obtained  so  far,  we  have  not  noticed  any  evidence  of 
numerically  introduced  discontinuity  of  the  solution  at  subdomain  interfaces.  Briley’s  outflow 
boundary  conditions  neither  introduce  artificial  boundary  layers  in  the  neighbourhood  of  the 
outflow  boundary  nor  do  they  affect  the  fiowfield  in  the  neighbourhood  of  the  cavity.  A  comment 
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relevant  to  the  wide  cavity  simulations  also  is  that  the  quality  of  the  solution  is  not  affected  in 
the  context  of  the  present  spectral  simulations  when  the  aspect  ratio  of  the  minimal  discrete  cell 
increases.  This  is  a  problem  which  is  of  concern  to  finite-difference  and  finite- volume  simulations, 
where  the  ratio  of  the  minimal  discretised  lengths  in  the  two  spatial  directions  must  be  of  0(1). 
Here  this  aspect  ratio  is  shown  to  take  values  from  0.25  to  10  inside  the  cavity  without  any 
visible  infiuence  on  the  quality  of  the  solution.  Furthermore,  the  pattern  predicted  by  Burggraf 
(1966)  on  the  basis  of  the  lid-driven  cavity,  including  its  symmetries,  is  also  present  in  our  deep 
open-cavity  numerical  results;  note  that  11  orders  of  magnitude  separate  the  maximum  and 
the  minimum  values  at  which  isolines  of  (  in  the  deepest  cavity  are  drawn  in  Figures  36-39 
and  40.  From  a  physical  point  of  view  we  observe  in  the  open  cavity  the  formation  of  multiple 
recirculation  zones  of  essentially  Stokes  fiow  (known  as  Moffat  eddies)  as  the  depth  of  the  cavity 
increases.  The  strength  of  recirculation  decreases  exponentially  as  one  progresses  from  y  =  0  to 
y  =  — 1£)|  with  most  of  the  activity  taking  place  in  the  neighbourhood  of  the  interface  between 
the  cavity  and  the  outer  flow.  From  the  point  of  view  of  a  subsequent  global  linear  instability 
analysis,  owing  to  their  small  magnitude  in  comparison  with  the  flow  in  the  neighbourhood  of 
the  cavity  lips,  the  details  of  the  Moffat  eddies  are  not  expected  to  influence  the  accuracy  of  the 
instability  results  to  be  obtained.  Nevertheless,  it  is  remarkable  that  a  resolution  of  the  flow  with 
20  collocation  points  per  spatial  direction  inside  the  cavity  is  capable  of  delivering  this  result. 
If  finer  resolution  of  the  Moffat  eddies  is  required,  this  may  be  straightforwardly  obtained  by 
refining  the  cavity  resolution  by  subdividing  it  in  additional  subdomains,  as  indicated  in  Figure 
3.  While  it  is  certainly  possible  to  pursue  higher  resolution  simulations,  we  bear  in  mind  the 
ultimate  objective  of  the  present  work,  which  is  recovery  of  as  accurate  as  possible  a  basic 
flow  on  as  small  number  as  possible  a  number  of  discretisation  points.  The  proposed  spectral 
multidomain  algorithm  has  been  shown  to  satisfy  this  requirement. 

(iii)  Solutions  in  a  full  bay  model 

Beyond  delivering  accurate  results  at  modest  resolutions,  the  multidomain  algorithm  proposed 
for  the  solution  of  the  open  cavity  problem  is  well-suited  to  solve  this  problem  beyond  the 
idealised  situation  of  an  empty  rectangular  open  cavity.  As  an  aside  to  the  main  theme  of  the 
present  report,  we  discuss  the  solution  to  a  problem  representative  of  the  class  of  flows  indicated 
in  §  2  (a )  in  which  an  object  is  contained  in  the  open  cavity.  A  two-dimensional  model  of  the 
domain  of  a  full  bay  is  shown  in  Figure  41.  This  domain  is  chosen  as  a  demonstrator  of  the 
capabilities  of  the  present  spectral  multidomain  algorithm,  without  any  pretension  to  model  a 
specific  cargo;  it  should  be  clear,  though,  that  the  limit  of  complexity  of  the  geometry  which  may 
be  tackled  using  our  approach  is  set  by  efficiency  aspects  and  the  number  of  processors  available 
to  be  devoted  to  the  (clearly  parallelisable)  solution  approach.  All  solutions  which  follow  were 
obtained  on  a  single-processor  platform  using  the  24  subdomains  indicated  in  Figure  41. 

From  an  algorithmic  point  of  view,  there  is  no  conceptual  addition  in  complexity  by  considering 
the  domain  shown  in  Figure  41.  The  nature  of  the  boundary  conditions  necessary  to  close  the 
problem  in  each  subdomain  is  the  same  as  that  in  the  empty  open  cavity,  namely  an  inflow 
boundary,  a  solid  wall  and  a  free-stream  or  outflow  boundary.  The  compatibility  conditions 
across  domain  interfaces  are  treated  as  part  of  the  new  iterative  algorithm  used  for  the  open 
cavity  solutions  presented  in  the  previous  sections.  The  efficiency  of  our  approach  is  indicated  by 
the  fact  that  solutions  on  8  x  8, 12  x  12, 16  x  16  or  20  x  20  conformable  grids  in  each  subdomain 
could  be  obtained  on  single-processor  platforms  within  reasonable  computing  times.  A  technical 
detail  of  interest  is  that  if  the  number  of  discretisation  points  is  kept  the  same  in  all  subdomains 
which  discretise  space  in  this  relatively  complex  geometry,  the  smallest  grid  spacing  in  the 
smallest  in  size  subdomain  imposes  an  unnecessarily  small  iteration  parameter  cj  (or  time-step  in 
a  time-accurate  framework)  if  standard  Chebyshev  Gauss-Lobatto  (CGL)  grids  are  used.  Indeed, 
the  value  of  u  scales  with  1/iV^,  where  N  is  the  number  of  CGL  points  used  in  each  spatial 
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direction;  this  limitation  severely  restricts  the  usefulness  of  a  spectral  method  as  resolution 
increases.  A  remedy  is  to  redistribute  the  CGL  points  in  a  more  equidistant  manner  compared 
with  the  (strongly  stretched  at  the  endpoints  of  the  integration  domain)  CGL  grid.  There  are 
several  ways  to  achieve  this  target;  here  we  have  constructed  a  straightforward  modification 
of  the  grid  proposed  by  Kosloff  and  TaLEzer  (1993).  Specifically,  we  map  the  CGL  points 
^  =  cos  =  0,  •  •  ■ ,  AT  onto 


,  XR-XL  ( arcsin[(^  -  1)^] 
^  2  1  arcsin(l  —  fl) 


(4.6) 


where  xl  and  xr  are  the  endpoints  of  the  original  domain,  mapped  on  themselves.  Appro¬ 
priate  choices  of  the  0[N~‘^)  parameter  fl  alleviate  the  l/N^  restriction  for  u)  to  l/AT;  using 
(jj  =:  10”^  we  have  obtained  converged  results  within  O(IO^)  iterations  in  what  follows.  The 
convergence  history  of  a  typical  solution  obtained  with  8  collocation  points  per  spatial  direction 
and  subdomain  is  shown  in  Figure  42.  The  large  number  of  iterations  shown  is  characteristic 
of  the  complexity  of  the  problem  solved.  While  the  discrepancy  between  successive  iterates  be¬ 
comes  lower  than  0.02%  after  ^  3  x  10^  iterations,  the  solution  relaxes  to  the  converged  result 
rather  slowly.  We  attribute  this  finding  to  the  relatively  large  fraction  of  subdomain  boundaries 
on  which  an  iteration  is  required  over  those  on  which  physical  boundary  conditions  are  supplied. 
This  fraction  increases  from  0.3  in  the  case  of  the  subdomain  discretisation  used  for  the  empty 
open  cavity  shown  in  Figure  2  (upper)  to  0.8  in  the  spatial  discretisation  of  Figure  41. 

Prom  a  physical  point  of  view,  the  objective  of  this  short  departure  from  the  main  theme  of 
our  report  has  been  to  compare  solutions  in  the  open  (empty)  and  the  present  model  of  a  full-bay 
cavity.  To  this  end,  we  chose  the  parameters  shown  in  Table  7  which  result  in  solutions  to  be 
compared  (in  the  sense  of  size  of  the  domain)  with  the  L/D  =  5  rectangular  cavity  results  also 
presented  in  Figure  34.  The  streamfunction  and  its  first  derivatives  are  shown  in  Figures  43-45 
along  the  results  of  the  equivalent-sized  empty  open  cavity.  In  all  results  an  equidistant  number  of 
contours  between  the  respective  maxima/minima  and  zero  are  provided.  As  expected,  substantial 
differences  in  the  solutions  in  the  two  configurations  exist;  the  presence  of  a  sizable  object  inside 
the  cavity  leads  to  the  destruction  of  the  zero-streamline  which  connects  the  upstream  and 
downstream  lips  of  the  empty  cavity  and  to  the  formation  of  strong  shear  layers  inside  the 
cavity,  with  new  stagnation  points  being  formed  as  a  consequence  of  the  chosen  geometry.  The 
extent  to  which  fluid  from  the  free-stream  entrains  into  the  open  cavity  may  be  seen  in  Figure 
43  where  the  regions  of  recirculating  fluid  under  the  object  have  been  highlighted. 

The  qualitative  features  of  the  results  presented  in  Figures  43-45  have  been  found  to  be  inde¬ 
pendent  of  resolution  when  using  between  8  and  20  points  per  spatial  direction  and  subdomain. 
Clearly,  running  the  code  on  a  parallel  machine  can  refine  the  quantitative  features  of  the  solu¬ 
tions  obtained;  this  is  also  beyond  the  scope  of  the  present  work.  What  we  are  interested  in  here 
is  viewing  these  results  within  the  framework  of  a  global  linear  instability  analysis.  The  conclu¬ 
sion  which  may  be  drawn  is  that  the  total  resolution  requirements  take  a  global  linear  analysis 
of  the  full-bay  basic  flow  model  well  beyond  the  reach  of  an  eigenvalue  problem  approach.  While 
it  might  still  be  possible  to  obtain  reasonably  well  resolved  global  linear  instability  results  in 
the  empty  open  cavity  using  an  eigenvalue  problem  approach  in  which  the  spatial  discretisation 
is  treated  in  a  coupled  manner,  a  time-accurate  approach  and  decoupled  calculation  of  spatial 
derivatives  is  essential  for  the  global  linear  instability  analysis  of  a  full-bay  model.  The  well- 
validated  spectral  multidomain  algorithm  proposed  herein  is  one  efficient  potential  candidate  to 
be  employed  in  this  respect. 
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==  Po'^o^T’o  =  have  been  used  for  the  nondimensionalisation  of  the  equations; 

these  have  been  combined  into  the  Reynolds  Re  =  po^io-^/Mo  and  Prandtl  Pr  =  numbers 

of  the  flow  with  the  latter  taken  to  assume  the  constant  value  for  air,  Pr  =  0.72. 

Stokes’  hypothesis 

A  =  (5.2) 

and  Sutherland’s  law  for  viscosity 


n  = 


1  +  Cs 

T  +  Cs' 


(5.3) 


with  C5  =  110.4°K/To  complete  the  system  of  equations  to  be  solved  in  the  framework  of 
compressible  direct  numerical  simulations. 


(6)  Inflow  boundary  conditions 

In  a  manner  analogous  to  the  incompressible  simulations  presented  in  §  4  at  the  inflow  bound¬ 
ary  in  the  DNS  the  similarity  profiles  pertaining  to  compressible  boundary-layer  flow  over  a 
flat  plate  at  given  free-stream  Reynolds  and  Mach  numbers  is  imposed.  The  similarity  solution 
is  obtained  using  the  classic  Howarth-Dorodnitsyn  transformation  of  the  compressible  laminar 
boundary  layer  equations  which  transforms  the  wall-normal  coordinate  y  into  a  similarity  coor¬ 
dinate  7]  defined  by 


where 


(5.4) 


V  = 


v/2^ 


and 


y 


(5.5) 


The  system  of  ordinary  differential  equations  to  be  solved  for  the  determination  of  the  stream- 
wise  velocity  component  u  and  temperature  T  is 


(x/")'+//"  =  0, 
(Xr')  +  fT'  +  {J-  i)Mlx  (/")'  =  0, 


(5.6) 

(5.7) 


where  f  —  u/uoo^X  ~  (pm)/(PooMoo)  and  primes  denote  differentiation  with  respect  to  the 
similarity  variable  rj.  The  system  (5.6-5. 7)  is  solved  subject  to  the  boundary  conditions 


/(O)  =  0, 

/(O)  =  0,  /(??  ->■  oo)  =  0, 

(5.8) 

alongside 

T{ri  -)•  oo)  =  1, 

(5.9) 

in  the  case  of  an  isothermal  wall,  or 
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t'{0)  =  0  (5.11) 

if  an  adiabatic  wall  is  considered.  Details  of  the  shooting  algorithm  used  and  solutions  of  the 
compressible  flat-plate  boundary  layer  obtained  are  outlined  in  the  Appendix. 

Solutions  of  the  compressible  flat-plate  boundary  layer  equations  have  been  obtained  using 
Pr  =  0.72,To  =  288.89^K  and  Moo  —  0.01,0.5,1.0  and  1.5.  Plots  of  the  streamwise  velocity 
component  u,  temperature  T  and  their  first  two  derivatives  are  shown  in  Figure  46.  Both  the 
viscous  and  the  thermal  boundary  layers  in  the  subsonic  cases  are  barely  distinguishable  from 
the  incompressible  result.  As  a  consequence,  a  (non-rational)  global  linear  instability  analysis 
performed  on  the  basis  of  the  incompressible  model  at  subsonic  Mach  numbers  and  the  consistent 
compressible  global  linear  instability  approach  are  expected  to  deliver  qualitatively  analogous 
results.  Of  interest  is  whether  the  presence  of  inflectional  profiles  in  the  compressible  basic  state 
at  inflow  introduces  additional  (global)  linear  instability  modes,  as  the  case  is  in  the  classic 
linear  analyses.  This  forms  one  of  the  questions  to  be  answered  in  future. 

(c)  Numerical  Methods 

After  discretising  the  open  cavity  domain  in  the  manner  presented  in  §1,  we  employ  in  each 
subdomain  a  twofold  extension  of  the  techniques  discussed  by  Theofilis  (1998b)  for  the  solution 
of  the  two-dimensional  compressible  Euler  equations.  Firstly,  in  the  present  problem  two  inhomo¬ 
geneous  spatial  directions  are  considered,  as  opposed  to  the  one  inhomogeneous  and  one  periodic 
spatial  direction  of  the  latter  work.  Secondly,  the  algorithm  is  straightforwardly  extended  to  in¬ 
corporate  calculation  of  the  viscous  fluxes  in  the  form  presented  above.  Spatial  derivatives  are 
calculated  using  the  same  collocation  algorithm  discussed  for  the  incompressible  simulations.  Al¬ 
ternative  schemes  based  on  compact  finite-difference  methods,  which  have  reached  a  high  level 
of  sophistication  (Visbal  and  Gaitonde  1998;  Gaitonde  and  Visbal  1999),  might  also  have  been 
used;  however,  the  issue  of  performance  of  the  numerical  scheme  for  the  calculation  of  spatial 
derivatives  was  found  to  be  highly  platform-dependent  (Theofilis  1998b)  and,  hence,  we  adhere 
to  the  spectral  scheme  which  provides  optimal  accuracy  at  modest  resolutions. 

By  contrast  to  the  incompressible  simulations,  here  we  concentrate  on  obtaining  time-accurate 
compressible  results  using  a  low-storage  fully  explicit  algorithm  due  to  Wray  (1986)  for  the  time- 
integration  of  the  equations  of  motion.  The  algorithm  advances  the  solution  vector  Q  according 
to 


ki  =  At  F(Q^t^) 

k2  =  At  F(Q^  -h  tn  +  2^^) 

ka  =  At  F(Q”  +  "ki  -h  Y^k2,  t^  -h  g^^) 

qn+l  ^  Qn  ^  1^1  +  ^k2  +  ^k3,  (5.12) 

where  superscript  (n)  indicates  time-level,  At  is  the  time-step  and  F  =  — Fc  —  Gc  +  l/ReF^  + 
1/ReGv,  as  defined  in  (5.1).  The  attractive  feature  of  the  particular  RK3  scheme  is  that  it 
provides  third-order  accuracy  in  time,  while  it  requires  only  two  levels  of  storage.  As  such  the 
RK3  represents  a  compromise  between  the  more  accurate  but  more  expensive  classic  RK4  and 
the  less  accurate  and  equally  expensive  second-order  accurate  RK2  scheme.  The  severe  time- step 
restrictions  of  explicit  schemes  in  conjuction  with  spectral  spatial  discretisation  are  alleviated  by 
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use  of  the  mapping  proposed  by  Kosloff  and  Tal-Ezer  (1993)  discussed  in  §  4  (d).  Appropriate 
choice  of  the  single  mapping  parameter  results  in  near-uniform  distribution  of  collocation  points 
in  the  two  spatial  directions,  which  in  turn  permits  time-steps  of  0{l/N)^  N  denoting  the  total 
number  of  collocation  points  in  either  the  x—  or  the  y— direction,  as  opposed  to  0(1/^^^)  by 
which  standard  Chebyshev  methods  are  restricted  (Canuto  et  al.  1987). 

A  straightforward  Dirichlet-Neumann  algorithm  has  been  constructed  to  ensure  continuity 
of  the  solution  at  an  interface  T  of  two  subdomains  fli  and  ^2,  according  to  which  a  solution  Qi 
of  (5.1)  is  obtained  in  fix  which  satisfies  an  arbitrary  Dirichlet  boundary  condition  and  provides 
Neumann  data  for  the  solution  Q2  of  the  governing  equations  in  $12-  The  algorithm  is 


(  _  ^Gi 

dt  dx  dy 

Q(m)  ^  ^(m-1) 

< 

dQ2  _  _  dG2 

dt  dx  dy 

,  dQ^^^/dnr  =  do!^^ /dnr 


in  Oi, 

on  r, 

in  ^2, 
on  r, 


(5.13) 


where  F  =  Fc  —  ;^F^;  and  G  =  Gc  —  The  vector  of  Dirichlet  data  on  the  interface  is 

updated  according  to 


aW  =  +  (1  -  w)  ,  (5.14) 

with  superscript  (m)  indicating  iteration  level  and  to  being  an  0(1)  relaxation  parameter. 
Note  that  the  iteration  produces  O^— continuous  values  of  the  conservative  quantities  at  each 
new  time-level,  from  which  the  primitive  variables  may  be  extracted. 


(d)  Boundary  conditions  for  aeroacoustics  calculations 

The  objective  currently  is  calculation  of  steady-state  solutions  in  an  open  cavity  in  compress¬ 
ible  flow.  Aside  from  the  issue  of  accuracy,  which  may  be  guaranteed  by  ensuring  sufficiently 
high  resolution,  a  novelty  introduced  by  compressibility  is  the  finite  speed  of  propagation  of 
pressure  waves  and  the  potential  of  confusion  between  acoustic  signals  generated  by  the  flow 
itself  and  numerically  generated  reflections  at  the  boundaries  of  the  calculation  domains.  This 
issue  is  clearly  demonstrated  by  the  problem  of  advection  of  a  vortical  disturbance  (Visbal  and 
Gaitonde  1998) 


p  =  1, 

(5.15) 

U  =  Uoo  -  -^{y  -  Vc)  exp  (-r^/2)  , 

(5.16) 

V  =  T(a;  _  Xc)  exp  , 

(5.17) 

(  2\ 

P  =  Poo  ^  j  ’ 

(5.18) 

where  (a;c,  yc)  is  the  location  around  which  vorticity  is  distributed  initially,  Uqq  is  the  advection 
velocity,  poo  is  the  ambient  pressure,  =  {x  —  Xc)‘^/R^-{-{y  —  yc)‘^/E?  and  C  and  R  are  constants. 
Two  sets  of  results  were  obtained  using  the  compressible  DNS  techniques  discussed;  the  results 
serve  to  illustrate  the  point  regarding  numerical  reflections  at  finite  boundaries.  All  simulations 
were  performed  in  single  domains  using  upwards  of  50  collocation  points  per  spatial  direction 
and  time-steps  of  0(10"^). 
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Firstly,  we  have  chosen  {xcUc)  =  (0,0),  [/qo  =  0,(7  =  l,i?  =  1  and  performed  several  sim¬ 
ulations  initialised  by  (5.15)  and  (5.18)  alone,  setting  u  =  v  =  0  initially;  the  transient  so¬ 
lution  quickly  settles  to  u  and  v  predicted  by  (5.16-5.17)  as  seen  in  Figure  47,  independently 
of  Reynolds  or  Mach  number  values.  While  no  distortions  are  to  be  seen  on  account  of  nu¬ 
merical  dissipation,  which  is  practically  absent  in  the  spectral  scheme  used,  these  transient 
solutions  are  destroyed  by  numerical  instability  as  time  in  the  simulation  progresses;  iden¬ 
tification  of  the  origin  of  this  numerical  instability  is  our  next  concern.  To  this  end  we  set 
(^ciVc)  =  (— 5,0),C/oo  =  1,(7  ==  0.1, R  =  0.75  and  perform  simulations  at  M  =  0.4, 0.8  and 
M  —  1.2  at  Re  =  10^,  initialising  the  calculation  using  (5.15-5.18).  The  results  at  the  high¬ 
est  Mach  number  are  presented  in  Figures  48  and  49.  If  attention  is  focussed  on  the  velocity 
components  alone  (or  vorticity,  not  presented  here),  one  sees  that  the  vortical  structure  is  being 
advected  by  the  fiow  at  the  imposed  speed  Uoo  throughout  the  simulation.  This  result  is  repeated 
at  all  Mach  numbers  examined,  the  difference  being  that  when  all  parameters  in  the  respective 
simulations  are  kept  identical  as  the  Mach  number  decreases  numerical  instability  sets  in  earlier 
in  time.  The  same  conclusion  is  to  be  drawn  from  the  results  of  Figure  50  where  the  passage 
of  the  vortex  is  recorded  at  two  fixed  locations  (x  =  a;o,y  =  0)  with  xq  =  0  and  xq  «  10, 
respectively.  Starting  at  Xc  =  —5  the  vortex  passes  from  these  locations  with  no  distortion,  as 
seen  clearly  in  the  u— signal  which  is  indistinguishable  from  the  analytical  solution  at  all  times 
examined. 

However,  the  reason  for  the  numerical  instability  is  clearly  to  be  seen  in  the  signals  of  p,  u 
and  V  as  well  as  in  the  density  and  pressure  isosurfaces  presented  in  Figures  48  and  49.  A  cylin¬ 
drical  pressure  wave  emanating  from  the  vortical  structure  propagates  in  space,  reflects  on  the 
boundary  and  returns  to  the  integration  domain.  The  reason  is  clearly  the  boundary  conditions 
imposed  at  the  artificially  truncated  boundaries,  the  reflecting  character  of  which  is  responsible 
for  the  numerical  instability  initially  and  the  destruction  of  the  simulation  eventually.  One  rem¬ 
edy  is  to  place  the  boundaries  far  enough  from  the  region  of  interest;  however,  this  approach  is 
not  suitable  for  simulations  in  the  open  cavity,  where  long-time  integrations  will  be  necessary  for 
a  steady  state  solution  to  be  obtained.  A  more  elegant  solution  to  this  problem  is  the  imposition 
of  non-reflecting  boundary  conditions  based  on  characteristic  analyses  (Thompson  1987;  Poinsot 
and  Lele  1992).  While  the  issue  of  non-reflecting  boundary  conditions  is  the  subject  of  current 
investigations  (Rowley  and  Colonius  2000)  it  should  be  noted  that  in  subsonic  simulations  some 
amount  of  reflection  must  be  allowed  in  order  for  pressure  to  be  able  to  adjust  to  the  ambient 
value  and  implementation  of  perfectly  non-reflecting  boundary  conditions  is  not  advisable.  Our 
current  objective  is  the  study  of  non-reflecting  far-field  boundary  conditions  in  conjunction  with 
spectral  multidomain  algorithms.  Specifically,  we  are  interested  in  the  interaction  of  pressure 
waves  with  artificially  created  internal  boundaries,  such  as  the  interfaces  of  the  multidomain 
algorithm,  and  the  potential  generation  of  instabilities  of  numerical  origin  at  such  boundaries. 
Once  this  issue  has  been  answered  in  a  satisfactory  manner  we  intend  to  proceed  and  obtain  the 
compressible  analoga  of  the  open  cavity  solutions  presented  in  §  4  (d);  results  will  be  presented 
in  due  course. 


Contract  No.  F61775-99-WE090 


Spectral  multidomain  for  laminar  flows  in  2D  open  cavities  25 

6.  Conclusions 

A  spectral  multidomain  algorithm  of  the  Dirichlet-Neumann  class  has  been  presented  for  the 
numerical  solution  of  the  Poisson  equation  and  of  the  streamfunction/vorticity- transport  formm 
lation  of  the  system  of  equations  governing  incompressible  fluid  flow,  when  space  is  decomposable 
in  rectangular  subdomains.  The  novelty  of  the  algorithm  is  that,  in  the  case  solutions  to  the 
steady  equations  of  motion  are  sought,  the  iteration  necessary  for  satisfaction  of  solution  continu¬ 
ity  across  subdomain  interfaces  and  that  for  the  nonlinearity  of  the  governing  partial  differential 
equations  are  combined  and  performed  in  a  single  step.  This  results  in  order-of-magnitude  sav¬ 
ings  compared  with  the  standard  approach  in  which  these  iterations  are  nested  within  each 
other.  The  issue  of  physical  boundary  conditions  and  numerical  compatibility  conditions  has 
been  discussed  by  reference  to  several  benchmark  solutions  with  which  excellent  agreement  has 
been  obtained.  In  geometries  in  which  single-domain  calculations  may  also  be  used  it  has  been 
observed  that  higher  resolution  is  necessary  when  using  multidomain  in  order  to  achieve  results 
of  the  same  quality  as  those  of  the  single-domain  calculation.  Relatively  high  resolutions  are 
necessary  for  the  solution  of  both  the  Poisson  and  the  incompressible  Navier-Stokes  and  con¬ 
tinuity  equations  when  the  domain  contains  geometric  singularities,  as  the  case  is  in  the  open 
cavity.  Particularly  revealing  are  the  solutions  of  the  Poisson  equation  in  such  domains,  where 
estimates  of  the  number  of  collocation  points  necessary  in  each  subdomain  in  order  for  solu¬ 
tions  to  converge  within  prescribed  tolerances  may  be  obtained.  This  is  especially  helpful  in  the 
context  of  the  global  linear  instability  analyses,  which  form  the  objective  of  an  extension  of  the 
present  work,  since  the  ability  to  increase  the  resolution  of  the  steady  laminar  basic  flow  at  will 
is  absent  in  the  framework  of  the  instability  analyses. 

Steady  laminar  incompressible  solutions  for  low  Reynolds  number  flow  in  arbitrary  aspect- 
ratio  open  cavities  have  been  obtained.  The  deep  cavity  limit  at  the  Reynolds  numbers  examined 
was  found  to  bear  remarkable  analogies  with  the  well-known  patterns  of  the  lid-driven  cavity 
flow  (Burggraf  1966).  Wide  open  cavity  solutions,  on  the  other  hand,  were  found  to  correspond 
to  the  union  of  the  well-known  flow  patterns  in  the  backward-  and  forward-facing  step  geome¬ 
tries.  The  high  resolution  requirements  suggested  by  the  model  problems  examined  earlier  were 
confirmed  in  the  open  cavity  simulations.  One  of  the  most  interesting  findings  of  the  present 
work  is  the  resonance-like  behaviour  of  the  flow  at  modest  spatial  resolutions.  A  clearly  defined 
periodic  solution  pattern,  which  may  well  be  confused  with  flow  resonance  between  the  upstream 
and  the  downstream  vertical  cavity  walls,  appeared  as  a  consequence  of  modest  resolutions  and 
disappeared  when  resolution  was  increased  further.  The  need  for  careful  further  numerical  ex¬ 
perimentation  in  order  for  the  critical  Reynolds  number  iie2d,st  for  onset  of  two-dimensional 
unsteadiness  has  thus  been  underlined.  A  further  interesting  result  of  the  present  work  is  the 
assessment  of  the  differences  of  the  flowfield  set  up  in  the  empty  as  opposed  to  an  open  cavity 
containing  an  object.  From  the  point  of  view  of  a  subsequent  global  linear  instability  analysis, 
the  resolution  requirements  in  latter  flow  were  found  to  be  such  that  the  instability  analysis 
may  only  be  performed  in  the  context  of  an  initial-boundary  value  problem.  The  resolution  re¬ 
quirements  of  the  empty  open  cavity,  on  the  other  hand,  are  such  that  a  partial-derivative  linear 
eigenvalue  problem  approach  (Theofilis  1998a)  may  be  employed. 

In  order  for  compressible  flow  solutions  in  the  open  cavity  to  be  obtained  the  issue  of  ap¬ 
propriate  boundary  and  compatibility  conditions  must  be  addressed.  The  problem  at  the  inflow 
boundary  was  solved  by  provision  of  the  compressible  flat-plate  boundary-layer  solution,  while 
the  potential  pitfalls  of  naive  application  of  the  incompressible  boundary  conditions  to  the  solu¬ 
tion  of  the  boundary  closure  problem  in  the  far-field  and  the  downstream  outflow  boundaries  of 
the  computational  domain  in  compressible  simulations  have  been  highlighted.  It  was  shown  that 
boundary  conditions  which  do  not  prevent  reflections  at  artificially  introduced  boundaries  may 
not  necessarily  destroy  the  simulation  but  can  definitely  lead  to  misinterpretation  of  its  results. 
We  are  currently  working  on  this  issue. 
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Our  recent  work  on  residuals  in  DNS,  motivated  by  our  desire  to  understand  the  origins  of 
instability  and  three-dimensionalisation  of  two-dimensional  steady  flows,  such  as  those  in  the 
open  cavity,  suggests  that  the  existence  of  two-dimensional  steady-state  solutions  in  the  open 
cavity  is  synonymous  with  stability  of  all  two-dimensional  global  flow  eigenmodes.  This  does 
not  prevent  three-dimensional  modes  being  unstable.  The  extent  to  which  the  conclusions  put 
forward  in  §  4  (d )  carry  weight  in  a  three-dimensional  environment  will  precisely  be  determined 
by  the  global  instability  analysis  in  which  the  entire  steady  (a:,  y)— flowfields  calculated  in  the 
framework  of  the  present  work  may  be  used  as  variable  coefficients  for  the  partial-differential- 
equation-based  complex  nonsymmetric  generalised  eigenvalue  problem.  Estimates,  based  on  the 
known  global  linear  instability  results  in  the  lid-driven  cavity,  of  the  instability  results  in  the 
open  cavity  may  already  be  obtained  indirectly  by  comparing  the  open  and  the  lid-driven  cavity 
basic  states;  this  work  is  also  in  progress. 

The  global  linear  instability  analysis  of  the  open  cavity  flow  forms  the  ultimate  objective  of  an 
extension  of  the  present  work;  an  intermediate  step  is  the  extension  of  the  incompressible  part 
of  the  present  work  to  implement  a  time-accurate  approach  for  the  recovery  of  open-cavity  solu¬ 
tions,  an  approach  which  is  expected  to  be  an  order  of  magnitude  more  intensive  computationally 
compared  with  the  iterative  algorithm  proposed  herein,  on  account  of  the  iterations  for  satis¬ 
faction  of  solution  continuity  across  subdomains  being  nested  within  the  fractional  time-steps 
of  the  time-integration  procedure;  however,  the  time-accurate  algorithm  is  necessary  in  order 
for  two-dimensional  laminar  time-periodic  basic  states  to  be  obtained.  Both  incompressible  and 
compressible  time-periodic  fields  may  form  an  alternative  basic  flow  compared  with  the  laminar 
steady-states  obtained  herein.  Qualitatively  diff’erent  basic  flows  and  a  global  linear  instability 
analysis  approach  incorporating  Floquet  theory  in  the  case  of  a  time-periodic  basic  state  are 
expected  to  result  in  the  identification  of  different  types  of  global  linear  instability  mechanisms. 
The  deliverable  of  this  approach  will  be  a  critical  Reynolds  number  i?e2d,tp  >  ^e2d,st  which 
will  be  closer  to  that  relevant  to  Air  Force  needs  in  comparison  with  i?e2d,st-  Clearly,  linear 
amplification  of  both  types  of  global  disturbances,  stationary  and  time-periodic^  is  relevant  to 
flight  Reynolds  numbers. 

Finally,  the  experience  obtained  in  the  present  work  suggests  that  the  resolution  necessary 
for  adequate  description  of  the  steady  laminar  basic  flow  in  the  open  cavity  at  high  Reynolds 
numbers  or  the  resolution  of  the  flowfield  set  up  by  the  presence  of  sizeable  objects  within 
the  cavity  calls  for  an  extension  of  the  partial-derivative-eigenvalue-problem-based  global  linear 
instability  approach  (Theofilis  1998a).  In  this  respect,  one  of  the  present  findings  points  at 
the  fact  that  a  solution  approach  in  which  all  subdomains  are  solved  in  a  coupled  manner 
is  impractical  on  account  of  the  large  total  number  of  collocation  points  and  the  associated 
memory  limitations.  One  way  forward  in  this  respect  is  implementation  of  an  initial-boundary- 
value  problem  approach  for  global  linear  instability  analysis  which  incorporates  the  validated 
spectral  multidomain  algorithm  discussed  in  the  present  work. 


Contract  No.  F61 775-99- WE090 


27 


Spectral  multidomain  for  laminar  flows  in  2D  open  cavities 

Appendix  A. 

A  shooting  algorithm  for  the  compressible  flat-plate  boundary  layer  equations 

Key  to  the  solution  of  the  system  (5.6-5.7)  is  the  boundary-layer  assumption  of  a  constant 
pressure  impressed  upon  the  flow  across  the  layer, 


It  follows  that  pT  =  1  and 


(Al) 


X^T2 


1  +  Cs 


(A  2) 


T  +  Cs’ 

where  Sutherland’s  viscosity  law  (5.3)  has  been  used.  Introduction  of  (A  2)  into  (5. 6-5. 7) 
permits  solution  of  this  system  for  the  similarity  variables  f{ri)  and  T{ri)  subject  to  boundary 
conditions  (5.8-5.9)  and  either  of  (5.10)  or  (5.11).  A  straightforward  shooting  approach  combined 
with  Newton  iteration  until  all  boundary  conditions  are  satisfied  is  employed  to  solve  the  system 


on  a  uniform  77— grid.  The  governing  equations  may  be  written  in  the  form  of  as  a  system  of 
ordinary  differential  equations 


where 


s', 


9i 

=  f 

,  hi 

=  92, 

92 

=  /' 

,  h2 

=  93, 

/ 

93 

=  /" 

,  hs 

=  -\9395 

9i 

=  T 

,  /14 

=  95, 

/ 

95 

=  t' 

,  h^ 

=  -\9l- 

:9i93, 


=  -^ff5-T9i55-Pr(7-l)Af^53- 


(A3) 


(A4) 


The  known  boundary  conditions  are  supplemented  with  estimates  of  the  unknown  values 
/(->  oo),/''(0),/"(^  oo),T'(->-  00)  and  either  of  T(0)  or  t'(O),  depending  on  whether  an 
isothermal  or  adiabatic  problem  is  considered,  respectively.  The  77-grid  encompasses  a  large 
number  of  uniformly-distributed  points,  such  that  subsequent  interpolations  of  the  similarity 
results  onto  the  DNS  grid  may  be  performed  without  appreciable  error  being  introduced. 
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Table  1.  Validation  of  the  multidomain  algorithm  for  the  Poisson  equation  in  the  rectangular  duct  steady  laminar 

flow;  x{y)  =  . 


(Nx.Ny) 

iterations 

CPU  time  (sec) 

'(pixL.O) 

8x8 

8 

0.14  (0.07) 

0.901343 

16  X  16 

6 

2.04  (0.73) 

0.901344 

32  X  32 

4 

73.21  (4.74) 

0.901344 

{Nx,Ny) 

iterations 

CPU  time  (sec) 

V’{XL,0) 

8x8 

8 

0.14  (0.07) 

0.955314 

16  X  16 

6 

2.04  (0.73) 

0.955321 

32  X  32 

4 

73.18  (4.99) 

0.955321 

^(0,0) 

1p{XR^0) 

IpxiXR,  0) 

0.152104 

0.981459 

0.967024 

-0.428399(-l) 

0.152106 

0.981459 

0.967035 

-0.428125(-1) 

0.152106 

0.981459 

0.967035 

-0.428125^1) 

‘lpx{XL,0) 

A  =  A 

mO) 

'IpiXR.O) 

1px{XR,0) 

0.699181  (-1) 

0.996145 

0.990301 

-0.140611(-1) 

0.699102(-1) 

0.996145 

0.990328 

-0.139337^1) 

0.699102(-1) 

0.996145 

0.990328 

-0.139337(-1) 
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8x8 

16  X  16 

32  X  32 

64  X  64 

^>(-1,0.5) 

V>a,(-1,0.5) 

^,,(-1,0.5) 

2.6(-15) 

0.817566 

7.9(-14) 

-1.8(-13) 

0.816394 

-3.4(-ll) 

-1.2(-11) 

0.816215 

7.6(-10) 

7.0(-12) 

0.816187 

-7.5(-9) 

V’(-0.75,0.5) 

V-i  (-0.75, 0.5) 
(-0.75, 0.5) 

0.153385 

0.452766 

-0.028354 

0.153105 

0.451631 

-0.027761 

0.153059 

0.451444 

-0.027659 

0.153052 

0.451413 

-0.027643 

^>(-0.5, 0.5) 
^>^,(-0.5, 0.5) 
i/>y(-0.5,0.5) 

0.242893 

0.282003 

-0.116120 

0.242337 

0.281842 

-0.101810 

0.242249 

0.281784 

-0.097664 

0.242235 

0.281769 

-0.096486 

^(0,0.5) 

1/11(0,0.5) 

V-y  (0,0.5) 

0.317491 

8.9(-6) 

-0.239739 

0.316856 

8.9(-7) 

-0.238431 

0.316751 

1.2(-7) 

-0.238214 

0.316735 

1.9(-8) 

-0.238179 

i/>(0.5,0.5) 

V'x(0.5,0.5) 

V'j,(0.5,0.5) 

0.242896 

-0.282063 

-0.076711 

0.242338 

-0.281845 

-0.090028 

0.242249 

-0.281784 

-0.094253 

0.242235 

-0.281769 

-0.095514 

V>(0.75,0.5) 

V-x  (0.75, 0.5) 
V>j,(0.75,0.5) 

0.153388 

-0.452788 

-0.028361 

0.153105 

-0.451633 

-0.027762 

0.153059 

-0.451444 

-0.027660 

0.153052 

-0.451413 

-0.027643 

i/.(l,0.5) 

V-*  (1,0.5) 
V-,(l,0.5) 

0 

-0.817611 

0 

0 

-0.816396 

0 

0 

-0.816216 

0 

0 

-0.816187 

0 
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Table  3.  Validation  of  the  multidomain  algorithm  for  the  Poisson  equation  in  the  backward-facing  step  geometry 

(Black  1997). 


16  X  16 

32  X  32 

64  X  64 

V>(0, 37r/4) 

-0.162463 

-0.161707 

-0.161585 

-0.279086 

-0.278571 

-0.278505 

^,(0.37r/4) 

-0.222365 

-0.234559 

-0.237928 

ip{n,0) 

0.645964 

0.645964 

0.645964 

V’a)(7r,0) 

1.2(-12) 

1.6(-11) 

1.5(-10) 

ipy(n,Q) 

-1.719610 

-1.719540 

-1.719530 
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Table  4.  Boundary  conditions  for  single- domain  calculations  in  the  PPF.  Dirichlet  and  Neumann  data  are 

denoted  by  (d)  and  (n)  respectively. 


Boundary 

Type 

c 

N 

(d) 

Eq.  (4.1) 

Eq.  (2.1) 

(n) 

1py=Q 

E 

(n) 

tpxx  —  0 

« 

II 

O 

S 

(d) 

Eq.  (4.1) 

Eq.  (2.1) 

(n) 

ipy=0 

w 

(d) 

Eq.  (4.1) 

Eq.  (2.1) 

(n) 

■(px  =0 
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Table  5.  Boundary  conditions  for  multidomain  calculations  in  the  square  lid- driven  cavity,  Dirichlet  and 
Neumann  data  are  denoted  by  (d)  and  (n)  respectively. 


Domain  1011 


Boundary 

xP 

c 

N 

II II 

o 

Eq.  (2.1) 

E 

(d) 

Eq.  (2.1) 

S 

i/;  =  0 
V’j/  =  0 

Eq.  (2.1) 

w 

i/;  =  0 

Eq.  (2.1) 

Domain  1110 


Boundary 

xP 

c 

N 

o 

II  II 

Eq.  (2.1) 

E 

II  II 

o  <= 

Eq.  (2.1) 

S 

^  =  0 
=  0 

Eq.  (2.1) 

w 

(n) 

Eq.  (2.1) 
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Table  6.  Convergence  of  open-cavity  solutions  at  different  Reynolds  numbers 
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Resolution 

50 

100 

Re 

500 

1000 

5000 

8x8 

0.091626 

0.090597 

0.079286 

0.072195 

0.139555 

12  X  12 

0.087852 

0.087006 

0.077489 

0.072145 

0.146231 

16  X  16 

0.085891 

0.085086 

0.076138 

0.071182 

0.147909 

20  X  20 

0.084899 

0.084103 

0.075375 

0.070629 

0.146547 
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Table  7.  Parameters  for  the  solution  of  the  fulhbay  open  cavity  model;  {xo^y^)  =  (0.5,0.121)  and  S*  is  the 

boundary  layer  thickness  at  inflow. 


{xi  -  xo)/S*  =  20 
{x2-xi)IS*  =  17 
(a;3  -  X2)/S"‘  =  6 
(x4  ~  xz)IS*  =  4 
{xs  —  X4)/S*  =  5 
(a:6  —  x^)l6*  —  10 
(xj  —  xq)I5*  =  13 
(xg  —  X'j^fd*  —  30 


2/4  =  0 

(2/4  -  yz)ld*  =  3 
(2/3  -  y2)IS*  =  2 

(y2  -  yi)ld*  =  2 
(2/1  -  2/o)/<^*  3 
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Table  8.  Convergence  of  the  transient  solutions  of  the  inviscid  vortical  flow  at  (Xyy)  =  (0,0)  and  t  =  0.5; 

w,z;<  10-^^ 


o 

t-H 

II 

,M  =  1.0 

iZe=  10^ 

,M  =  0.6 

ile  =  10^ 

II 

p 

E 

P 

E 

P 

E 

25  X  25 

0.605211 

3.16587 

0.633437 

5.26767 

0.547016 

2.40012 

50  X  50 

0.614062 

3.19235 

0.636089 

5.28667 

0.560949 

2.34479 

100  X  100 

0.614119 

3.19272 

0.636078 

5.28658 

0.561028 

2.34295 
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Figure  1.  The  building  blocks  of  the  two-dimensional  domains  considered 
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Figure  2.  Spatial  decomposition  of  an  empty  (upper)  and  an  open  cavity  containing  a  protrusion  (lower)  using 

the  building  blocks  of  Figure  1. 
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Figure  3.  Example  of  spatial  resolution  refinement  by  decomposing  an  original  (upper)  into  further  subdomains 

(lower) . 
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Figure  4.  The  streamfunction  in  flow  in  an  aspect-ratio  A  =  3  rectangular  duct  (Tatsumi  &  Yoshimura,  199( 
recovered  using  three  subdomains  with  interfaces  at  xl  =  — 0.5A  and  xr  =  0.25A. 
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Figure  5.  The  streamfunction  ?/?  in  a  grooved  channel  geometry  recovered  using  five  subdomains 
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Figure  6.  The  velocity  component  u  in  a  grooved  channel  geometry  recovered  using  five  subdomains 
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Figure  7.  The  velocity  component  v  in  &  grooved  channel  geometry  recovered  using  live  subdomains 


Contract  No.  F61 775-99- WE090 


Spectral  multidomain  for  laminar  flows  in  2D  open  cavities 


45 


-0.8 


-0.6  -0,4  -0.2  0 


X 


•e  8.  Solution  of  =  1  subject  to  the  boundary  C( 


99-WE090 


onditions  of  Black  (1997). 


46 


V.  Theofilis 


Vjfo.o)  C|,(0,0) 


Figure  9.  Spectral  multidomain  solution  of  the  plane  channel  laminar  basic  flow.  Shown  are  the  flow  streamfunction 
(left),  the  streamwise  velocity  component  u  (centre)  and  the  wall- normal  derivative  of  the  flow  vorticity  C  (right) 
functions  of  the  number  of  iterations  N  and  the  relaxation  parameter  r  at  (x^y)  —  (0,0).  r  =  0(0.01)  and 
r  =  0(0.1)  was  used  during  the  early  and  late  iteration  cycles,  respectively. 
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Figure  10.  The  streamwise  velocity  component  and  the  vorticity  in  plane  Poiseuille  flow  solved  using  two 

domains  whose  interface  is  at  x  =  0. 
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Figure  11.  The  streamfunction  and  vorticity  (upper)  and  the  velocity  components  (lower)  in  the  square  lid-driven 
cavity  at  Re  —  100.  Solution  obtained  using  48  collocation  points  per  spatial  direction  in  each  of  two  domains 
whose  interface  is  at  a;  =  0.5. 
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Figure  12,  Upper:  Convergence  history  of  ip  (left)  and  C  (right)  as  function  of  the  time^like  variable  r.  Lower: 
Collapse  of  ip  (left)  and  ipy  (right)  on  the  respective  Blasius  profiles.  The  imposed  at  inflow  Blasius  solutions  are 
indicated  by  the  solid  lines;  results  obtained  at  the  interface  of  the  domains  are  shown  by  (o)  and  those  at  the 
outflow  boundary  are  denoted  by  (□). 
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Figure  13.  Instability  of  the  numerical  scheme  in  Blasius  flow  at  low  resolutions  and  high  w  values. 
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Figure  14.  The  streamfunction  and  vorticity  (upper)  and  the  velocity  components  (lower)  in  incompressible  flow 
over  a  flat-plate.  Solution  obtained  using  16  collocation  points  per  spatial  direction  in  each  of  the  two  subdomains 
connected  at  rCc  =  0.7. 
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Figure  16.  The  streamfunction  and  vorticity  (upper)  and  the  velocity  components  (lower)  in  incompressible  flow 
over  a  flat-plate.  Solution  obtained  using  48  collocation  points  per  spatial  direction  in  each  of  the  two  subdomains 
connected  at  Xc  =  0.21. 
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Figure  17.  Collapse  of  (left)  and  ipy  (right)  on  the  respective  Blasius  profiles.  The  imposed  at  inflow  Blasius 
solutions  are  indicated  by  the  solid  lines;  results  obtained  at  the  interface  of  the  domains  are  shown  by  (o)  • 
and  those  at  the  outflow  boundary  are  denoted  by  (□).  Upper:  Calculation  using  16  x  16  collocation  points  per 
subdomain.  Lower:  Calculation  using  48  x  48  collocation  points  per  subdoman. 
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Figure  18.  Schematic  representation  of  the  geometry  and  definition  of  parameters  in  the  open-cavity  flow. 
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Figure  20.  Steady**state  solutions  at  Re  =  100  obtained  with  20^  collocation  points  in  each  subdomain. 
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Figure  21.  Steady-state  solutions  at  Re  =  500  obtained  with  20^  collocation  points  in  each  subdomain 
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N(o 


Figure  25.  Oscillatory  open  cavity  solutions  obtained  on  account  of  poor  spatial  resolution  using  16^  colloca¬ 
tion  points  in  each  subdomain.  Shown  are  the  values  of  ipx^'ipy  and  ^  at  the  midpoint  of  the  interface  between 
subdomains  II  and  IV. 


Contract  No,  F61775-99-WE090 


Figure  26.  Upper:  transient  oscillatory  solutions  of  the  velocity  components  in  the  open  cavity  corresponding  to 
the  result  of  Figure  25.  Lower:  steady-state  solutions  of  the  same  quantities  obtained  with  20^  collocation  points 
in  each  subdomain. 
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Figure  34.  Steady-state  rectangular  open  cavity  solutions.  Left  to  right  column:  Upper  to  lower  row: 

LjD  =  1,356,7,9  with  D  =  10(5*  and  6*  the  flat-plate  boundary  layer  displacement  thickness  at  the  inflow 
boundary.  A  spatial  discretisation  comprising  20^  collocation  points  in  each  subdomain  was  used. 
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Figure  40.  Steady-state  rectangular  open  cavity  solutions.  Left  to  right  column:  V'l/j  V'®:  C*  Upper  to  lower  row: 
D/L  =  1, 2, 3, 4  with  D  =  and  <5*  the  flat-plate  boundary  layer  displacement  thickness  at  the  inflow  boundary. 
A  spatial  discretisation  comprising  20^  collocation  points  in  each  subdomain  was  used. 
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Figure  41.  A  model  of  an  object  contained  in  the 
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Figure  42.  Local  iteration  history  of  the  maxima  of  C  and  the  velocity  components  for  flow  over  the  configuration 
depicted  in  Figure  41  using  the  parameters  of  Table  7.  Re  =  295  is  imposed  at  inflow  and  20^  collocation  points 
are  used  in  each  of  the  24  subdomain  shown  in  Figure  41. 
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Figure  44.  Velocity  component  ^ 
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Figure  46.  Compressible  flat-plate  boundary  layer  similarity  solutions.  Upper:  Streamwise  velocity  component 
u  —  f  and  its  first  two  derivatives;  Lower:  Temperature  T  distribution.  Solid:  Moo  =  0.01,  dotted:  Moo  =  0.5, 
dashed:  Moo  =  1.0,  dash-dotted:  Moo  =  15. 
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Figure  49.  Advected  vortex  in  =  10^, M  =  1.2  viscous  compressible  flow  (Visbal  and  Gaitonde  1998).  Left  to 
right  column  p^u^v^p;  upper  to  lower  row,  t  —  12.5, 15.0, 17.5 
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